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ABSTRACT 

Currently  there  is  no  generally  accepted  procedure  for  calculation  of  structural 
loadpaths,  which  wotild  show  how  remote  loads  are  equilibrated  through  a  structure 
and  could  provide  insight  into  how  well  a  structure  is  performing  its  intended  load- 
canying  functions.  Kelly  and  Elsley  have  recently  proposed  a  method  for  computing 
loadflow  orientations  and  loadpaths  using  finite  element  results,  which  is  based  on 
iterative  solutions  of  non-linear  equations.  In  this  paper,  we  have  enhanced  their 
theoretical  formulation  and  general  procedure  by  deriving  explicit  expressions  for 
computing  loadflow  orientations.  The  new  equations  produce  more  accmate  loadflow 
orientations  compared  to  the  prior  approach  and  improve  the  fidelity  of  calcxflated 
loadpaths.  In  a  series  of  benchmark  problems,  we  have  investigated  non-optimal  and 
optimal  holes  in  plates  using  loadflow  visualisation  to  identify  their  key  features.  We 
found  that  recirctdation  is  apparent  for  non-optimal  hole  shapes,  whereas  no 
recirculation  zone  is  present  for  optimal  shapes.  Although  very  highly  refined  finite 
element  meshes  were  utilised,  the  implications  are  that  even  more  refined  meshes  are 
required  to  fully  capture  the  complex  behaviour  that  exists  in  recirculation  zones.  The 
removal  of  the  recirculation  zone  for  a  non-optimal  shape  leads  to  a  better  shape,  but 
the  improvement  in  peak  stress  is  insignificant.  The  calculation  of  loadflow 
orientations  using  the  new  equations  is  simple,  and  could  be  used  with  any  firute 
element  analysis  code,  while  a  plotting  package  is  required  to  display  loadpaths. 
Loadflow  visualisation  is  a  powerful  tool  for  use  by  structural  designers  to  improve 
their  understanding  of  structural  performance,  the  application  of  which  can  potentially 
result  in  worthwhile  improvements  in  structural  efficiency. 
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Executive  Summary 

The  use  of  loadflow  visualisation  to  assist  the  structural  analyst  in  establishing  the 
performance  of  structures  under  load  is  an  emerging  new  field.  Recent  developments 
in  methodology  and  computer  implementation  of  loadflow  visualisation,  through  the 
use  of  finite  element  airalysis  (FEA)  results,  offer  the  potential  for  useful  engineering 
insights  into  the  way  a  given  structure  is  performing  its  intended  function.  This  holds 
promise  for  the  use  of  loadflow  visualisation  techniques  in  the  evaluation  of  airframe 
performance,  with  expected  applications  to  life  assessment  and  the  development  of  life 
extension  strategies  for  aircraft.  Hence,  loadflow  visualisation  may  be  particularly 
valuable  for  highly  stressed  airframes,  leading  to  longer-term  outcomes  such  as 
increased  residual  life  and  cost  benefits  due  to  longer  inspection  intervals. 

Currentiy  there  is  no  generally  accepted  procedure  for  calculation  of  structural 
loadpaths,  which  wotild  show  how  remote  loads  are  equilibrated  through  a  structure 
and  could  provide  insight  into  how  well  a  structure  is  performing  its  intended  load- 
carrying  hmctioris.  KeUy  and  Elsley  have  recently  proposed  a  method  for  computing 
loadflow  orientations  and  loadpaths  using  finite  element  results,  which  is  based  on 
iterative  solutions  of  non-hnear  equations.  In  this  paper,  we  have  enhanced  their 
theoretical  formulation  and  general  procedure  by  deriving  explicit  expressions  for 
computing  loadflow  orientations.  The  new  equations  produce  more  accurate  loadflow 
orientations  compared  to  the  prior  approach  and  improve  the  fidelity  of  calculated 
loadpaths.  In  a  series  of  benchmark  problems,  we  have  investigated  non-optimal  and 
optimal  holes  in  plates  using  loadflow  visualisation  to  identify  tiieir  key  features.  We 
found  that  recirculation  is  apparent  for  non-optimal  hole  shapes,  whereas  no 
recirculation  zone  is  present  for  optimal  shapes.  Although  very  highly  refined  finite 
element  meshes  were  utilised,  the  implications  are  that  even  more  refined  meshes  are 
required  to  fully  captuxe  the  complex  behaviour  that  exists  in  recirculation  zones.  The 
removal  of  the  recirculation  zone  for  a  non-optimal  shape  leads  to  a  better  shape,  but 
the  improvement  in  peak  stress  is  insignificant.  The  calculation  of  loadflow 
orientations  using  the  new  equations  is  simple,  and  could  be  used  with  any  finite 
element  analysis  code,  while  a  plotting  package  is  required  to  display  loadpaths. 
Loadflow  visualisation  is  a  powerful  tool  for  use  by  structural  designers  to  improve 
their  understanding  of  structural  performance,  the  application  of  which  can  potentially 
result  in  worthwhile  improvements  in  structural  efficiency. 

The  developments  in  loadflow  visualisation  reported  here  provide  structural  analysts 
in  Airframes  and  Engines  Division  (AED)  with  the  capability  and  the  tools  leading  to  a 
better  imderstanding  of  the  performance  of  structural  components.  This  will  enable 
more  comprehensive  analysis  of  critical  components  in  RAAF  aircraft,  resulting  in 
potential  extensions  to  fatigue  life  and  increased  inspection  intervals. 
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1.  Introduction 

The  use  of  loadflow  visualisation  to  assist  the  structural  analyst  in  establishing  the 
performance  of  structures  under  load  is  an  emerging  new  field.  Loadflow  visualisation 
can  be  considered  to  consist  of  two  distinct  categories:  i)  display  of  loadflow 
orientations  at  discrete  points  in  a  body,  and  ii)  display  of  loadpaths  (trajectories) 
taken  through  a  body  by  unit  quantities  of  load.  Recent  developments  in  methodology 
and  computer  implementation  of  loadflow  visualisation,  tiirough  the  use  of  finite 
element  analysis  (FEA)  results,  offer  the  potential  for  useful  engineering  insights  into 
the  way  a  given  structure  is  performing  its  intended  function.  This  holds  promise  for 
the  use  of  loadflow  visualisation  techniques  in  the  evaluation  of  airframe  performance, 
with  applications  to  life  assessment  and  the  development  of  life  extension  strategies 
for  aircraft,  such  as  shape  optimisation  of  reworks  [1].  Thus,  loadflow  visualisation 
may  be  particularly  valuable  for  highly  stressed  airframes,  leading  to  longer-term 
outcomes  such  as  increased  residual  life  and  cost  benefits  due  to  longer  inspection 
intervals. 

There  is  a  very  limited  amount  of  published  prior  work  on  loadflow  and  its 
visualisation.  In  the  past,  the  orientations  of  the  major  principal  stresses  throughout  a 
body  have  often  been  used  to  try  to  depict  the  flow  of  load  through  a  structure. 
Making  use  of  photo-elasticity  in  an  early  paper.  Baud  [2]  noted  that  the  mapping  of 
principal  stress  direction  lines  in  structures  tiiat  have  stress  concentrations  could  be  a 
substantial  aid  in  explaining  stress  concentration.  However,  it  has  been  clearly 
demonstrated  by  Kelly  and  Elsley  [3]  [4]  that,  when  used  to  represent  loadflow 
orientations,  the  principal  stress  vectors  can  have  an  angular  error  of  up  to  45°  in 
regions  of  high  shear  stress. 

In  other  work.  Chaperon  et  al.  [5]  have  used  a  displacement-based  formulation  for 
loadflow  visualisation  that  is  based  on  the  orientation  of  energy-flow  vectors  to  study 
optimal  shapes,  but  neither  loadpath  lines  nor  compaction  of  loadpaths  were 
investigated.  An  optimal  shape  is  defined  as  one  that  minimises  the  peak  stress 
concentration  factor  (SCF)  in  the  structure  for  the  chosen  loading  condition.  They 
identified  "recirculation"  zones  and  considered  that  these  could  be  a  feature  of  a  non- 
optimal  hole  shape  with  redundant  material  present  that  might  possibly  be  removed 
with  only  a  minimal  effect  on  the  loadflow  in  the  structure.  This  concept  may  assist  in 
producing  an  initial  design  for  subsequent  shape  optimisation  of  the  structure. 
However,  tiie  energy-flow  vector  approach  has  shown  itself  to  be  sensitive  to  the 
location  of  reference  axes,  and  it  has  problems  in  handling  S5munetry  conditions, 
which  are  routinely  utilised  to  reduce  the  computational  size  of  finite  element  (FE) 
models.  Follow-on  work  by  Ibrahim  et  al.  [6]  and  Jones  and  Chen  [7]  on  the  use  of 
energy  density  vectors  has  overcome  the  symmetry-related  problems  that  complicate 
the  use  of  energy-flow  vectors. 

Furthermore,  much  of  the  available  literature  is  predominantly  concerned  with 
qualitative  descriptions  of  loadflow  without  any  method  of  calculation.  In  two 
engineering  publications  dealing  extensively  with  stress  concentration  factors. 


1 


DSTO-RR-0166 


Peterson  [8]  and  Pilkey  [9]  refer  to  a  flow  analogy  when  discussing  aspects  of  stresses 
around  notches.  However,  the  pictures  of  "flow  lines"  in  and  around  notches, 
shoulder  fillets,  T-heads  and  nuts  are  purely  of  a  schematic  nature.  The  concept  of 
"force  lines,"  where  the  spacing  of  the  flow  lines  may  be  regarded  as  being  inversely 
proportional  to  the  intensity  of  the  stress,  is  similarly  utilised  by  Osgood  [10]  to  help 
visualise  the  effects  of  notches  and  fillets  on  stress  distribution.  Broek  [11]  also 
considers  loadpath  (loadflow)  lines  around  notches,  holes,  and  fillets  where  these 
discontinuities  cause  a  deviation  in  the  loadflow  lines,  leading  to  stress  concentration. 
In  his  work  on  minimising  stress  concentrations  in  fillets.  Baud  [12]  used  a 
hydrodynamical  fluid  flow  analogy  to  show  that  a  streamline  fillet,  whose  contour  is 
the  same  as  the  contotu  for  a  free  jet  of  water  issuing  from  a  circular  orifice,  produces  a 
low  SCF.  However,  Neuber,  in  his  treatment  of  the  theory  of  notch  stresses  [13],  has 
used  the  integration  of  stress  components  to  anal5dically  calculate  the  lines  of  stress 
arovmd  an  internal  notch  for  a  case  where  an  analytical  solution  was  feasible. 

Kelly  and  Elsley  [3]  [4]  have  recently  made  significant  inroads  into  the  visualisation  of 
loadflow.  They  have  developed  a  theoretical  basis  for  calculating  loadflow  orientations 
at  discrete  points  within  general  structures  by  iterative  solution  of  non-linear 
transcendental  equations.  Their  method  relies  on  using  the  results  from  a  FEA, 
computed  at  the  nodes  and/  or  element  centroids  in  the  FE  mesh.  Kelly  and  Elsley 
demonstrated  that  their  method  is  accurate  in  regions  of  high  shear  stress,  unlike 
loadpath  visualisation  based  on  the  use  of  principal  stress  trajectories.  Along  the  lines 
of  Neuber  [13],  they  implemented  an  automated  procedtue  for  integrating  the  stresses 
obtained  from  the  FEA  solution  for  determining  loadpath  contours,  and  they 
successfully  applied  this  to  a  cracked  plate.  In  later  work,  Kelly  [14]  demonstrated  how 
loadflow  visualisation  could  assist  in  defining  an  efficient  FE  mesh,  where  the 
orientation  of  the  mesh  is  aligned  with  the  loadpaths  through  the  structure.  Kelly  and 
Tosh  [15]  have  also  studied  the  potential  for  using  loadpath  lines  to  help  wiA  the 
optimal  placement  of  fibres  in  fibre-reinforced  composite  laminates.  In  recent 
tmpublished  work,  Kelly  and  Tosh  [16]  have  bypassed  the  inconvenience  of  the 
integration  method  of  computing  loadpaths  by  using  directly  the  vector-like  field  of 
loadflow  orientations  calculated  from  the  FEA  stress  values. 

Initially,  this  paper  reports  on  a  key  result  extending  Kelly  and  Elsley's  [3]  [4]  loadflow 
visualisation  procedtire  by  developing  explicit  expressions  for  calculating  loadflow 
orientation,  and  these  aie  given  in  Section  2.  The  derivation  presented  here  also 
disperses  a  degree  of  confusion  arising  from  reference  to  the  major  and  minor 
principal  stress  angles  in  the  presentation  of  the  original  theory.  In  Section  3,  the 
loadflow  orientations  and  loadpaths  arotmd  non-optimal  and  optimal  rework  shapes 
are  investigated  using  the  new  simplified  equations,  which  have  obviated  the  need  for 
an  iterative  numerical  solution.  This  is  an  area  that  has  not  previously  been  considered 
in  depth.  By  using  much  more  refined  FE  meshes  than  those  analysed  by  other 
workers,  it  was  possible  to  perform  loadflow  visualisation  at  a  greater  level  of  detail 
than  has  previously  been  the  case.  In  producing  the  loadpath  lines,  the  new  method 
proposed  by  Kelly  and  Tosh  [16]  was  applied,  and  a  detailed  description  of  the 
me&odology  is  given  here.  The  use  of  the  compaction  of  loadpath  lines  to  compute 
stress  concentration  factors  is  also  demonstrated.  In  Section  4  the  potential  for  using 
recirculation  zones  that  are  present  in  non-optimal  shapes  to  lead  to  first  estimates  of 
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Optimal  shapes  for  use  in  shape  optimisation  procedures  is  investigated.  A  discussion 
and  some  conclusions  are  presented  in  Section  5. 


2.  Theory 

2.1  Analysis  of  loadflow  orientations 

Kelly  and  Elsley  [3]  [4]  have  developed  a  procedure  that  can  be  applied  to  tiie 
calctdation  of  loadflow  orientations  within  a  structure,  using  the  results  obtained  from 
a  FEA.  We  will  present  the  theory  from  the  start  and  make  significant  improvements 
along  the  way.  Although  Kelly  and  Elsley's  logic  will  be  retained,  we  will  make  use  of 
the  notation  ti^iat  is  commonly  found  in  standard  texts  on  the  theory  of  elasticity  in 
order  to  make  the  derivation  presented  here  easier  to  understand. 

Kelly  and  Elsley  defined  a  loadpath  to  be  the  trajectory  taken  within  a  structure  by  a 
tmit  quantity  of  appUed  load,  beginning  at  a  point  of  application  and  ending  at  an 
equilibrating  reaction.  They  introduced  tire  concept  of  a  hypothetical  force  "stream 
tube",  which  is  depicted  in  Figure  1.  For  any  selected  resolved  force  direction,  depicted 
as  the  ;r-direction  in  Figure  1,  there  is  no  flow  of  load  across  the  boundary  of  the  force 
tube,  and  tiie  loadpath  is  therefore  bounded  by  lines  along  which  there  is  no 
contribution  to  the  force  in  the  ;r-direction.  Because  equilibrium  must  be  satisfied  over 
the  tube,  we  have  tiiat 

F,a-F,b  =  0  (1) 

Note  that  in  Figure  1  the  applied  and  reaction  loads  Fxa  and  Fxb  are  aligned  with  the  x- 
direction,  but  an  equivalent  diagram  can  also  be  produced  for  the  forces  in  the  y- 
direction  (or  any  other  chosen  direction). 

On  a  plane  that  is  tangential  to  the  force  tube  wall  there  are  two  components  of  stress 
acting,  a„  and  .  Here  is  the  local  stress  perpendicular  to  the  plane  and  is  the 

local  shear  stress  acting  along  the  plane,  where  is  oriented  at  an  angle  0  with 
respect  to  a  reference  x-axis,  as  depicted  in  Figure  1.  It  should  be  noted  that  the  angle 
G  defined  here  is  different  to  that  used  by  KeUy  and  Elsley,  and  the  derivation  will  be 
presented  from  the  start  for  clarity. 

The  condition  of  force  equilibriiim  is  satisfied  locally  along  any  segment  of  loadpath 
wall  if  the  force  components  due  to  the  action  of  c„  and  in  the  direction  parallel  to 
the  applied  loads  cancel  out.  For  appHed  forces  resolved  in  the  a:-direction,  the  value 
0  =  0^  that  satisfies  this  condition  is  given  by  the  solution  to  the  equation 

a„  cos  0,  -  sin  0^  =  0  (2) 
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For  applied  forces  resolved  in  the  y-direction,  the  eqtiivalent  equation  that  must  be 
satisfied  to  determine  the  value  0  =  0^  is 

cT„sin0^+T„,cos0^  =0  (3) 

From  the  standard  two-dimensional  theory  of  elasticity  [17],  the  following  equations 
define  the  normal  stress  a„  and  the  shearing  stress  at  a  point  on  an  arbitrary  plane 
oriented  at  an  angle  0  with  respect  to  a  reference  r-axis  in  terms  of  the  known  stresses 
(j^,  and  ,  as  depicted  in  Figure  2. 

a„  =  cos^  0  +  sin^  0  +  sin  0  cos  0  (4) 

=  -(cr^  ~  )  sin  0  cos  0  +  (cos^  0  -  sin^  0)  (5) 

In  passing,  it  can  be  noted  that,  by  substituting  0  =  0°  and  0  =  90°  into  Equations  (4) 
and  (5),  it  is  possible  to  obtain  the  original  stresses  cj^,  and  x^ .  These  act  on  the 

original  x-  and  y-planes,  which  are  orthogonal,  as  shown  in  Figure  2.  In  Kelly  and 
Elsle5^s  work.  Equations  (2)  and  (3)  were  treated  as  being  non-linear  and  were 
subsequently  solved  in  an  iterative  manner,  and  the  stresses  and  Xj^  were 

obtained  by  FEA.  However,  it  will  be  shown  that  recoiurse  to  a  non-linear  solution 
technique  is  unnecessary,  as  some  additional  algebraic  manipidation  presents  us  with 
explicit  analytical  expressions  for  the  loadflow  orientation  angles,  as  computed  in  the 
two  chosen  orthogonal  directions. 

Apart  from  straightforward  differences  in  nomenclature,  it  should  be  noted  tiiat  in  [3] 
and  [4]  the  expression  for  x„,  has  what  appears  to  be  a  typographical  error.  In 
Equation  (3)  therein,  the  term  involving  is  missing  a  minus  sign,  although  the 
loadflow  orientations  that  were  presented  in  those  two  references  are  correct. 

Substituting  the  expressions  for  and  x„,  in  Equations  (4)  and  (5)  into  Equation  (2), 
and  after  some  algebraic  manipulation  and  simplification,  we  finally  obtain 

tan0^=-^  (6) 

In  a  similar  manner.  Equation  (3)  becomes 

tan0  =-^  (7) 

It  can  be  seen  that  Equation  (6)  for  computing  the  r-direction  loadflow  orientations  is 
independent  of  o^,  and  Equation  (7)  for  computing  the  y-direction  loadflow 
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orientations  is  independent  of  .  This  characteristic  was  not  previously  apparent  in 
prior  work.  At  any  given  point  along  the  force  tube  wall,  the  angle  a  of  the  loadflow 
orientation  is  given  by  the  tangent  to  the  force  tube,  as  depicted  in  Figure  1.  The  angle 
a  can  therefore  be  expressed  in  terms  of  0  as 

a  =  0  +  90°  (8) 

Thus,  at  any  point  along  the  force  tube  wall,  the  local  loadflow  orientations  and 
tty  of  the  resolved  x-direction  and  y-direction  loadpath  lines  are  simply 

a,=0,+9O°  (9) 

a^=  0^+90°  (10) 

Using  the  above  two  relations,  the  equations  for  the  actual  loadflow  orientations 
and  Uy  become 


tana,=^  (11) 

tana  =-^  (12) 

Thus,  the  derivation  presented  here  has  led  to  the  development  of  explicit  equations 
for  the  x-direction  and  y-direction  loadflow  orientations,  a  result  that  has  not 
previously  been  available.  Equations  (11)  and  (12)  can  easily  be  solved  to  produce  a, 
and  tty  at  any  point  in  a  two-dimensional  stress  field.  This  key  result  considerably 

simplifies  the  computation  of  loadflow  orientations  throughout  a  body,  compared  to 
the  procedure  requiring  the  solution  of  non-linear  transcendental  equations  by  an 
iterative  numerical  procedaire  that  was  used  in  [3]  and  [4].  This  makes  it  particularly 
easy  to  compute  the  loadflow  orientations  from  the  stresses  produced  from  a  FEA  of  a 
structure. 

2.2  Principal  stress  directions 

For  the  purpose  of  comparing  loadflow  orientations  with  principal  stress  directions,  it 
is  convenient  to  include  here  the  standard  definition  of  the  well-known  major  and 
minor  principal  stress  angles,  0^i  and  9^2'  respectively,  which  correspond  to  the 

major  and  minor  principal  stresses,  and  0^3  •  The  two  principal  stress  angles  are 
determined  using  the  two  solutions  to  the  equation 
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2t 

tan2e  = - ^ 

a. -o 


(13) 


Note  that  the  two  angles  0^,  and  0^2  are  90°  apart,  where  one  value  of  0^  always  lies 
between  ±45°  (inclusive),  and  the  other  one  is  90°  greater.  The  two  principal  stresses 
can  be  computed  by  substituting  0  =  0pj  and  0  =  0^2  “^to  Equation  (4),  and  the  major 

principal  stress  tr^j  corresponds  to  the  algebraically  larger  of  the  two  stresses. 


2.3  FEA  method  for  displaying  loadflow  orientations  and  loadpaths 

Although  no  theory  has  been  given  at  this  juncture,  Kelly  and  Tosh  [16]  have  put 
forward  a  method  for  taking  a  set  of  computed  loadflow  orientations  and 
subsequently  using  the  Tecplot  [18]  commercial  data  visualisation  program  to  plot 
loadpaths.  They  have  proposed  that  loadpaths  can  be  considered  to  be  analogous  to 
streamtraces,  where  a  streamtrace  is  the  path  traced  by  a  massless  particle  placed  at  an 
arbitrary  location  in  a  steady-state  vector  field. 

The  method  used  here  takes  the  nodal  stresses  computed  from  a  standard  FEA 
solution,  together  with  the  element  topologies,  and  stores  these  data  in  files  for 
subsequent  use.  Appendix  A  gives  a  det^ed  description  of  the  method  that  was  used 
when  working  with  the  MSC.NASTRAN  FEA  program  and  MSC.PATRAN  pre-  and 
post-processor,  together  with  our  custom-written  Fortran  utility  programs  used  to 
create  and  store  the  required  data.  These  are  used  to  extract  nodal  and  element 
information  from  the  NASTRAN  data  deck,  as  well  as  in  combining  the  nodal  and 
element  stress  data  extracted  using  PATRAN.  The  full-field  stress  results  are  then 
processed  using  our  own  in-house  custom-written  Fortran  program,  which  is 
described  and  listed  in  Appendix  B,  and  the  loadflow  orientations  at  each  node  in  the 
original  FE  mesh  are  computed  using  the  equations  developed  in  Section  2.1.  (The 
principal  stress  directions  are  also  calculated  for  comparison  purposes.)  The  results  are 
then  stored  as  a  field  of  unit  vectors  in  a  Tecplot-compatible  format  and  are 
subsequently  read  in  by  the  Tecplot  program.  Using  the  concept  tiiat  loadflow  lines  are 
analogous  to  streamtraces  [16],  Tecplot' s  built-in  facility  for  plotting  streamtraces  was 
utilised  to  compute  and  display  loadpaths  for  tihe  structures  of  interest.  The  detailed 
step-by-step  procedure  is  presented  m  Appendix  C. 


3.  Benchmark  Problems  For  Holes  in  Plates 

The  FE  models  of  square  plates  with  holes  that  were  analysed  here  were  created  and 
meshed  in  MSC.PATRAN.  The  MSC.NASTRAN  FEA  program  was  then  used  to 
perform  a  linear  elastic  static  analysis  to  compute  the  stresses  in  the  structures  of 
interest,  although  other  general  purpose  FEA  computer  software  could  also  have  been 
used.  The  computer  that  was  used  was  a  dual-processor  Hewlett-Packard  K260  server. 
Four-noded  quadrilateral  finite  elements  were  used  for  all  analyses,  which  were 
conducted  imder  two-dimensional  plane  stress  conditions.  The  material  properties 
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were  those  for  alitminium,  with  Young's  modulus  E  =  70x103  MPa  and  Poisson's  ratio 
V  =  0.30.  In  all  the  FE  models  analysed  here,  quarter  symmetry  was  used  for  the 
displacement  restraints,  which  consisted  of  Mx  =  0  along  the  y-axis  and  My  =  0  along  the 
x-axis.  Although  only  one  quarter  of  the  plate  could  have  been  modelled  in  the  FEA 
(due  to  symmetry),  it  was  decided  to  use  complete  plate  models  in  order  to  retain 
flexibility  in  the  post-processing  and  displa)dng  of  results.  One  of  the  aims  was  to 
allow  a  full  view  of  the  loadpaths  through  Ihe  structure  if  desired. 

The  meshes  that  were  used  were  very  refined,  particularly  near  the  holes,  in  order  to 
accurately  calculate  the  stresses  in  areas  of  high  stress  gradient.  In  contrast,  the  FE 
meshes  used  by  Kelly  and  Elsley  [3]  [4]  were  at  least  an  order  of  magnitude  coarser.  In 
any  case,  our  FE  models  are  not  aH  that  large  by  modem  standards,  and  therefore  their 
solution  poses  no  computational  difficulties.  By  using  highly  refined  FE  meshes  it  was 
possible  to  resolve  the  fine  details  in  recirculation  zones  that  would  have  been  missed 
had  coarser  meshes  been  utilised. 

3.1  A  circular  hole  in  a  uniaxial  stress  field 

3.1.1  Geometry  and  loading 

The  plate  geometry  and  loading  arrangement  that  was  considered  for  this  analysis  is 
shown  in  Figxire  3.  Here  the  200  mm  square  plate  contains  a  centrally  located  circular 
hole  with  20  mm  diameter.  A  remote  uniaxial  100  MPa  tensile  stress  field  aligned  in 
the  y-direction  is  applied  to  the  plate  along  the  plate  edges  at  y  =  ±100  mm. 

This  case  is  considered  to  be  a  good  test  problem  because  the  maximum  stress 
concentration  factor  at  the  hole  will  be  approximately  3.0,  and  therefore  the  hole  is 
expected  to  significantly  disturb  the  loadpath  lines  passing  in  its  vicinity.  The  size  of 
the  plate  is  considered  to  be  large  in  comparison  to  the  dimensions  of  the  hole,  so  the 
effect  of  the  finite  plate  width  is  expected  to  be  small.  The  peak  stress  concentration 
factor  will  occur  at  two  points  aroimd  the  hole,  located  at  {x,  y)  =  (±10,  0).  Another 
feature  of  this  geometry  is  that  there  is  a  region  of  compressive  botmdary  stress  at  the 
top  and  bottom  of  the  hole,  and  the  magnitude  of  this  compressive  stress  is  similar  in 
value  to  fiiat  of  the  applied  load.  The  complete  mesh  discretisation  is  shown  in  Figure 
4a,  and  consists  of  9200  elements.  Figure  4b  shows  a  close-up  view  of  the  mesh  arotmd 
the  circular  hole,  where  the  highly  refined  nature  of  the  mesh  is  evident. 

3.1.2  Principal  stress  directions 

The  orientations  G^i  of  the  principal  stress  aroxmd  the  circular  hole  are  shown  in 
Figure  5a,  while  Figure  5b  shows  the  orientations  of  the  principal  stress  ^ 

these  two  figures  the  angles  of  orientation  of  the  principal  stresses  are  represented  by 
the  orientations  of  the  individual  short  line  segments.  Vertical  line  segments  signify 
principal  stress  angles  of  90“,  horizontal  line  segments  represent  angles  of  0°,  while 
line  segments  orientated  at  other  angles  represent  cases  lying  between  these  two 
conditions. 
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In  Figure  5a  it  is  evident  that  the  principal  stresses  are  largely  aligned  at  an  angle 
of  90°,  which  corresponds  to  the  direction  of  the  applied  loading.  The  orientation  of 
closely  follows  the  boundary  of  the  hole  but,  as  the  top  of  the  hole  is  approached 

in  an  anti-clockwise  direction,  the  orientations  twist  arotmd  until  they  are  largely 
perpendicular  to  the  boundary  of  the  hole.  A  similar  effect  is  noted  in  Figure  5b,  but 
this  time  the  orientation  of  the  principal  stress  closely  follows  the  hole  boundary 

at  the  top  of  the  hole,  but  it  becomes  aligned  perpendicvdar  to  the  hole  botmdary  as  we 
move  in  a  clockwise  sense  from  the  top  to  the  right-hand  side  of  the  hole.  In  both 
cases,  the  transition  points  for  these  changes  in  the  orientations  of  the  principal 
stresses  and  occur  at  the  same  point  on  the  boimdary  of  die  hole. 

It  should  be  noted  that,  depending  on  the  location  along  the  free  boimdary  of  the  hole, 
at  least  one  of  the  two  principal  stresses,  CTp,  and  /  will  be  identically  equal  to  zero 

(normal  to  the  boundary).  The  non-zero  principal  stress  will  be  aligned  such  that  it  is 
tangential  to  die  boundary.  At  the  four  locations  where  the  hoop  stress  goes  from 
being  tensile  in  nature  to  being  compressive,  they  will  both  be  zero.  Even  when  a  very 
refined  mesh  is  used  in  the  FE  model,  these  conditions  will  usually  only  be 
approximated,  due  to  the  inherent  numerical  error  (although  small)  of  the  FEA.  Hence, 
in  Figure  5a  and  Figure  5b,  regions  where  the  principal  stress  orientations  are  shown 
as  being  perpendicular  to  the  hole  boimdary  are  in  fact  regions  where  that  particular 
principal  stress  is  actually  zero. 

3.1.3  Loadflow  orientations 

In  direct  comparison  to  the  principal  stress  orientations.  Figure  6a  shows  the  resolved 
r-direction  loadflow  orientations,  and  Figure  6b  shows  the  resolved  y-direction 
loadflow  orientations.  Qose-up  views  of  the  resolved  r-direction  and  y-direction 
loadflow  orientations  around  the  circular  hole  are  shown  in  Figure  7a  and  Figure  7b.  In 
a  qualitative  sense,  the  y-direction  loadflow  orientations  in  Figure  6b  are  the 
counterpart  of  the  principal  stress  orientations  of  Figure  5a,  and  vice  versa.  Both  of 

the  loadflow  orientation  patterns  are  considerably  more  complex  than  those  for  the 
principal  stress  orientations.  In  Figure  7b  it  can  clearly  be  seen  that  near  the  top  of  the 
hole  there  is  a  distinct  "recirculation"  region  present,  to  use  a  fluid  flow  analogy, 
where  the  loadflow  orientations  quite  distinctly  vary  through  360°.  Looking  along  the 
hole  boimdary  in  Figure  7b,  it  is  evident  that  tiie  notional  boundary  of  the  recirculation 
zone  begins  at  approximately  the  same  position  as  the  transition  point  described  in 
relation  to  the  orientations  of  principal  stress  shown  in  Figure  5a. 

3.1.4  Loadpath  contours 

The  resolved  jc-direction  and  y-direction  loadpaths  around  the  circular  hole  are  shown 
in  Figure  8a  and  Figure  8b.  The  loadpaths  were  produced  using  the  streamtrace 
plotting  feature  available  in  Tecplot.  At  this  stage  there  is  no  direct  theoretical  basis  as 
to  why  this  technique  could  be  used.  However,  it  appears  that  there  is  some  sort  of 
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analogy  between  streamtraces  in  fluid  flow  and  ihe  loadpath  theory  that  is  presented 
here. 

The  resolved  ;c-direction  loadpaths  in  Figure  8a  consist  of  closed  cxirves.  From  this 
characteristic,  it  appears  that  the  a:-direction  loadpaths  are  self-equilibrating,  and  there 
is  no  net  flow  of  load  in  this  resolved  direction.  The  loadpaths  were  obtained  by 
defining  six  streamtrace  rakes  radiating  outwards  from  the  origin,  oriented  at  angles  of 
±90°,  ±22.5°,  and  ±157.5°  relative  to  the  ±-axis.  Except  for  the  rakes  at  ±90°,  which  were 
terminated  at  one  end  at  {x,  y)  =  (0,  ±100),  each  rake  was  fully  contained  within  the 
plate  perimeter,  i.e.  did  not  terminate  anywhere  on  the  plate  edge  or  hole  boundary . 
Five  streamtraces  were  initiated  from  each  rake.  A  minor  degree  of  "spiralling'  is 
evident  in  some  of  the  loadpaths  (see  the  darker  lines  in  the  loadpaths  aligned  at  157.5° 
and  -22.5°).  Loadpath  spiralling  occurs  whenever  a  loadpath  swirls  inwards  upon 
itself  with  ever  decreasing  radius,  rather  than  producing  a  closed  contour.  In  order  to 
reduce  the  effects  of  spiralling  for  this  plot,  it  was  necessary  to  decrease  the  maximum 
step  size  used  by  the  Tecplot  streamtrace  plotting  integration  algorithm  from  the 
default  value  of  0.25  of  the  element  size  to  0.125.  Also,  the  maximum  number  of 
integration  steps  was  reduced  from  ihe  default  value  of  1000  to  450  in  order  to  reduce 
the  total  pathlength  of  each  of  the  loadpaths. 

For  the  resolved  y-direction  loadpaths  in  Figure  8b,  the  start  and  end  points  of  the 
loadpath  lines  were  chosen  to  be  uniformly  distributed  along  the  top  and  bottom 
edges  of  the  large  square  plate.  This  results  in  equal  quantities  of  load  between 
loadpath  lines.  The  loadpaths  are  largely  vertically  aligned  over  most  of  the  plate,  but 
they  diverge  around  the  hole,  indicating  that  the  hole  acts  as  an  unpediment  to  the 
flow  of  load  between  the  top  and  bottom  of  the  plate.  Looking  at  the  loadpath  lines 
that  pass  close  to  the  hole,  it  can  be  seen  that  they  first  start  to  diverge  at  a  distance  of 
approximately  three  hole  diameters  from  the  hole  centre.  Looking  closely  at  the 
spacing  of  the  loadpath  lines  that  run  across  the  horizontal  axis  passing  tiirough  the 
centre  of  the  hole,  it  is  evident  that  there  is  a  small  deflection  of  the  loadpath  lines 
towards  the  hole  (this  has  also  been  depicted  in  [15]).  This  is  opposite  to  the  overall 
general  behaviour  of  the  loadpath  lines  that  spread  out  away  from  the  hole  in  order  to 
pass  around  it. 

Figure  9a  and  Figure  9b  show  enlarged  views  of  the  resolved  x-direction  and  y- 
direction  loadpaths,  respectively,  in  the  vicinity  of  the  recirculation  zone  that  is  present 
near  the  top  of  the  hole  (see  Figure  7b).  In  order  to  minimise  the  effects  of  spiralling 
when  producing  Figure  9a,  the  maximum  integration  step  size  was  set  to  0.15,  and  the 
maximum  number  of  steps  was  set  to  250.  Also,  each  of  the  loadpaths  was  started  just 
outside  the  visible  plot  area  for  additional  control.  Looking  at  the  two  loadpaths  near 
the  top  of  the  hole  in  Figure  9b,  spiralling  is  evident  for  Ihe  imterminated  contour  that 
runs  around  the  closed  contour.  After  considerable  manual  intervention,  this  was  the 
best  that  could  be  achieved  for  this  plot,  even  after  adjustment  of  Tecplot  streamtrace 
integration  step  size  and  total  number  of  steps.  Because  die  mesh  is  large  in  relative 
terms  compared  to  the  detail  that  needs  to  be  shown,  it  is  anticipated  that  additional 
mesh  refinement  in  this  area  would  assist  in  producing  better  loadpath  contours  in  this 
small  region  that  is  exhibiting  complex  loadflow  recirculation  behaviour.  It  has  also 
been  our  experience  that,  when  originally  using  the  non-linear  solution  method  to 
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obtain  the  loadflow  orientations,  the  lower  accuracy  adversely  affected  the  production 
of  the  loadpath  contours  in  some  regions  of  the  structure,  sometimes  causing  the 
loadpaths  to  be  prematurely  terminated. 

3.1.5  Calculation  of  Kt from  loadpath  compaction 

When  looking  at  the  loadpath  contours  in  Figure  8b,  it  is  apparent  that  regions  where 
the  loadpaths  are  compacted  together  are  indicative  of  areas  that  are  experiencing 
higher  stress  levels.  By  comparing  the  relative  spacing  (i.e.  compaction)  of  the 
loadpath  contours  in  the  vicinity  of  the  hole  to  the  spacing  at  the  edges  of  the  plate 
where  the  load  is  applied,  it  is  possible  to  calculate  the  SCF  at  different  distances  from 
the  hole.  This  enables  an  independent  check  of  the  loadpath  contours  to  be  conducted. 

Figure  10a  shows  a  circular  hole  of  diameter  2fl  in  a  square  plate  that  is  loaded  along 
the  top  and  bottom  edges  by  a  uniform  stress  S.  For  this  same  loading  condition. 
Figure  10b  shows  an  elemental  force  tube  whose  width  is  centred  on  a  location  x  = 
Xb  at  the  top  and  bottom  edges  of  the  plate.  As  it  passes  around  the  hole,  the  force  tube 
compacts  and  its  width  becomes  centred  on  a  location  x  =  Xc  along  tiie  x-axis.  The 
SCF  at  the  location  x  =  Xc  along  the  x-axis  is  given  by 

SCF(a:.)  =  4^  (14) 

Ac 

For  tile  uniaxially  loaded  square  plate  with  a  circular  hole,  and  using  Equation  (14), 
the  SCF  values  computed  from  loadflow  compaction  measurements  at  different 
positions  along  the  x-axis  from  the  edge  of  the  hole  to  the  edge  of  the  plate  are  given  in 
Table  1.  When  performing  these  calculations,  tiie  force  tube  width  along  the  x-axis  was 
chosen  first,  and  then  the  loadpath  lines  were  computed  and  their  width  measured  at 
tiie  loaded  edges  of  the  plate.  This  permitted  better  control  of  the  loadpath  hnes  to 
enable  them  to  be  calculated  very  close  to  the  edge  of  the  hole.  Table  2  gives  the  ^F 
values  that  were  computed  from  nodal  stresses  in  the  FE  model. 

Taking  these  two  sets  of  results.  Figure  11  compares  the  SCF  computed  from 
compaction  of  loadpath  hnes  with  that  obtained  from  nodal  stresses  in  the  FE  model.  It 
is  apparent  from  these  results  that  the  SCF  as  computed  from  loadpath  compaction  is 
in  very  good  agreement  with  that  obtained  from  the  nodal  stresses.  As  a  result  of  the 
finite  width  of  the  plate,  the  value  of  the  SCF  at  the  edge  of  the  hole  computed  from 
the  FE  results  is  3.12.  (It  should  be  noted  that  for  a  circular  hole  in  an  infinite  plate  the 
expected  theoretical  value  of  SCF  is  3.00.)  The  SCF  obtained  at  the  plate  edges  at  (x,  y) 
-  (±100, 0)  is  0.96,  also  as  a  result  of  finite-width  effects. 

3.2  A  circular  hole  in  a  4:1  biaxial  stress  field 

This  particrilar  problem,  consisting  of  a  circular  hole  in  a  y:x  =  4:1  biaxial  stress  field, 
was  used  because  the  circular  hole  is  clearly  a  non-optimal  shape  for  tiiis  loading 
condition,  and  we  were  interested  in  contrasting  the  results  with  those  for  an  optimal 
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hole  at  a  later  stage.  It  is  well  known  that  an  elliptical  shape  in  an  infinite  plate  is 
optimal  for  this  case,  significantly  reducing  the  SCF,  and  the  elliptical  hole  case  is 
subsequently  studied  in  Section  3.3. 

3.2.1  Geometry  and  loading 

The  plate  geometry  and  loading  arrangement  that  was  considered  for  this  analysis  is 
shown  in  Figure  12.  Here  the  200  mm  square  plate  contains  a  circular  centrally  located 
hole  of  20  mm  diameter.  A  remote  4:1  biaxial  stress  field  is  applied  to  the  plate,  where 
the  applied  tension  stresses  are  25  MPa  along  the  edges  x  =  ±100  mm  and  100  MPa 
along  the  edges  y  =  ±100  mm.  A  feature  of  this  geometry  is  tiiat  there  is  a  region  of 
compressive  boundary  stress  at  the  top  and  bottom  of  the  hole.  The  highly  refined 
mesh  shown  earlier  in  Figure  4  was  once  again  utilised  for  the  FEA.  For  an  infinite 
plate  the  peak  stress  at  the  edge  of  the  hole  at  (r,  y)  =  (±10, 0)  is  known  to  be  275  MPa, 
indicating  that  there  is  a  severe  stress  concentration  effect  caused  by  the  presence  of 
the  hole.  The  resrdts  from  the  FEA,  which  include  finite-width  effects,  gave  a  slightly 
greater  peak  stress  of  284.7  MPa  (3.5%  higher).  Also,  the  stress  around  the  boundary  of 
the  hole  is  markedly  non-uniform,  which  indicates  that  by  choosing  a  more  optimal 
shape  the  stress  concentration  effect  of  the  hole  could  be  reduced  [1]. 

3.2.2  Loadflow  orientations 

Figure  13a  shows  the  resolved  r-direction  loadflow  orientations,  where  one  small  and 
distinct  recirculation  zone  is  present  and  is  located  at  the  top  of  the  hole.  This  is  in 
contrast  to  the  uniaxial  case,  which  has  two  large  regions  of  recirculation  in  this 
quadrant  of  the  plate.  Figure  13b  shows  the  resolved  y-direction  loadflow  orientations, 
and  a  small  region  of  recirculation  can  be  identified.  This  is  in  contrast  to  the  uniaxial 
case  where  a  very  distinct  recirculation  zone  was  present.  Although  the  largest 
component  of  the  applied  load  is  that  which  is  applied  in  the  y-direction,  it  is 
noteworthy  that  it  is  the  picture  of  the  x-direction  loadflow  orientations  that  more 
distinctly  shows  the  presence  of  a  recirculation  zone. 

3.2.3  Loadpath  contours 

The  resolved  x-direction  loadpaths  around  the  circular  hole  are  shown  in  Figure  14a, 
while  the  resolved  y-direction  loadpaths  are  depicted  in  Figure  14b.  In  both  cases  the 
loadpaths  are  uniformly  distributed  along  opposing  edges  of  the  square  plate,  as 
appropriate  for  the  loading  orientation  being  considered. 

In  Figure  14b  the  y-direction  loadpaths  are  vertically  aligned  fliroughout  a  large 
proportion  of  flie  plate  in  a  manner  that  is  very  similar  to  what  occurred  for  the 
uniaxially  loaded  plate  (see  Figure  8b).  In  comparison  to  this,  the  x-direction  loadpaths 
in  Figure  14a  show  considerably  more  deviation  from  the  horizontal  over  a  large  part 
of  flxe  plate.  The  curvature  of  the  x-direction  loadpath  lines  arormd  the  hole  is  much 
more  pronounced  than  for  the  y-direction  loadpaths,  and  is  quite  pronounced  at  even 
at  three  or  four  hole  diameters  away  from  the  centre  of  the  hole.  Looking  at  the  x- 
direction  loadpaths  immediately  above  and  below  the  hole,  it  is  quite  noticeable  that 
they  are  angled  in  towards  the  horizontal  centreline  of  the  plate,  starting  right  at  the 
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left  and  right  edges  of  the  plate.  Figure  15  shows  a  close-up  view  of  the  loadpaths  in 
the  vicinity  of  the  recirculation  region. 

3.3  An  optimal  elliptical  hole  in  a  4:1  biaxial  stress  field 

This  particular  geometry,  consisting  of  an  elliptical  hole  in  a  y:jc  =  4:1  biaxial  stress 
field,  was  chosen  because  it  is  well  known  that  an  elliptical  hole  is  the  optimal  shape 
for  this  loading  condition  in  an  infinite  plate,  in  the  sense  that  it  achieves  the  lowest 
value  of  SCF.  The  stress  field  that  exists  aroimd  such  an  elliptical  hole  produces  an 
almost  uniform  hoop  stress  around  the  boxmdary,  and  the  existence  of  a  state  of 
uniform  stress  along  botmdary  segments  is  a  feature  of  optimal  shapes  [1]. 

3.3.1  Geometry  and  loading 

The  plate  geometry  and  loading  arrangement  that  was  considered  for  this  analysis  is 
shown  in  Figure  16.  Here  the  square  plate  contains  a  centrally  located  elliptical  hole 
with  its  major  axis  aligned  with  the  y-direction.  The  major  axis  of  the  elliptical  hole  is 
80  mm  in  length  and  its  minor  axis  is  20  mm  in  width.  The  dimensions  of  the  plate  and 
tire  remote  4:1  biaxial  stress  field  are  the  same  as  those  used  for  the  circular  hole  case  in 
the  previous  section.  The  complete  mesh  discretisation  for  this  problem  is  shown  in 
Figxire  17a.  TTie  mesh  is  highly  refined  around  the  boundary  of  the  hole,  and  a  close-up 
view  is  provided  in  Figure  17b. 

For  an  infiiute  plate,  the  theoretical  peak  stress  for  this  loading  is  125  MPa,  and  the 
stress  is  completely  tiniform  arotmd  the  hole  botmdary.  This  peak  stress  is  much  lower 
than  that  obtained  for  a  circular  hole  in  the  same  stress  field,  which  produces  a 
theoretical  peak  stress  of  275  MPa.  It  is  interesting  to  note  that  the  peak  stress  for  the 
elliptical  hole  is  considerably  lower  even  though  the  area  of  the  elliptical  hole  is  much 
greater  than  that  of  tiie  circular  hole  that  was  studied  previously.  The  computed 
stresses  around  the  boundary  of  the  elliptical  hole  from  the  FE A,  which  includes  finite- 
width  effects,  were  133.7  MPa  ±8.2%  (a  mean  stress  of  128.8  MPa  with  a  standard 
deviation  of  ±3.6%  of  the  mean).  The  greatest  variation  from  the  theoretical  solution 
occurs  in  regions  near  the  top  and  bottom  of  the  eUipse,  where  the  hole  boimdary 
comes  closest  to  the  plate  edges.  In  comparison,  the  stresses  on  the  left  and  right 
boxmdaries  of  the  ellipse  at  and  around  the  region  where  x  =  ±10  mm  are  within  ±1% 
of  the  theoretical  values  for  the  infinite-plate  case.  Although  the  results  indicate  that 
the  elliptical  hole  in  a  finite-width  plate  that  is  being  analysed  here  is  sub-optimal  to 
some  degree,  the  difference  between  it  and  the  true  optimal  shape  for  this  loading  and 
plate  geometry  are  not  considered  to  be  very  large.  Hence,  it  is  felt  that  the  present 
elliptical  hole  in  a  finite-width  plate  can  still  be  assumed  to  be  optimal  for  our 
purposes. 

3.3.2  Loadflow  orientations 

The  resolved  x-direction  loadflow  orientations  are  shown  in  Figure  18a,  while  the 
resolved  y-direction  loadflow  orientations  are  shown  in  Figure  18b.  In  both  of  these 
cases  the  loadflow  orientations  follow  the  local  hole  botmdary  very  closely.  Unlike  the 
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case  for  the  non-optimal  circular  hole,  there  is  no  identifiable  recircidation  zone  in 
eidier  of  the  plots  of  loadflow  orientations,  which  indicates  the  optimality  of  the 
elliptical  hole. 

3.3.3  Loadpath  contours 

Figure  19a  and  Figure  19b  show  the  resolved  r-direction  and  y-direction  loadpaths, 
respectively.  These  were  produced  by  uniformly  distributing  the  loadpath  starting 
points  along  tiie  plate  edges  corresponding  to  r  =  100  and  y  =  100.  Some  extra 
tiniformly  distributed  :c-direction  loadpaths  were  added  to  the  area  where  the 
loadpaths  were  flowing  in  the  vicinity  of  the  elliptical  hole.  In  botii  parts  of  Figure  19, 
the  loadpaths  closely  hug  the  boundary  of  the  elliptical  hole  over  a  great  portion  of  its 
length.  The  fact  that  the  elliptical  hole  perturbs  the  y-direction  loadpaths  less  than  does 
a  circular  hole  is  another  indication  of  tire  optimality  of  the  elliptical  hole.  Looking  at 
tihe  r-direction  loadpaths  in  Figure  19a,  it  is  apparent  that,  in  the  vicinity  of  the 
horizontal  z-axis,  they  diverge  away  from  the  horizontal  as  they  approach  the  centre  of 
the  ellipse.  In  contrast  to  this,  the  loadpaths  in  Fig.  14(a)  that  pass  close  to  the  circular 
hole  actually  converge  towards  the  horizontal  centreline  in  the  vicinity  of  the  hole. 


4.  Prediction  of  Initial  Rework  Shapes 

In  the  previous  section,  a  number  of  examples  have  been  presented  that  show  the 
existence  of  loadflow  recirculation  regions  for  non-optimal  rework  shapes,  while 
optimal  rework  shapes  do  not  have  a  recirculation  region.  This  leads  to  the  h3^othesis 
that  those  areas  where  recirculation  occurs  can  be  taken  to  be  indicators  of  subregions 
in  a  structure  tiiat  can  be  removed  to  achieve  a  more  optimal  shape.  In  view  of  tiiis,  it 
was  decided  to  investigate  the  possibility  of  using  the  limits  of  the  recirculation  zone 
to  define  the  boundary  of  an  initial  rework  region  that  may  provide  a  good  starting 
shape  for  the  application  of  a  shape  optimisation  procedure. 

4.1  Reworked  circular  hole  in  a  uniaxial  stress  field 

The  loadflow  orientations  obtained  around  a  circular  hole  in  a  large  square  plate 
uniaxially  loaded  in  the  y-direction  have  previously  been  shown  in  Figure  6  and 
Figure  7.  The  most  significant  recirculation  zone  corresponds  to  the  resolved  y- 
direction  loadflow  orientations  (Figure  7b),  and  the  rework  region  is  selected  to 
coincide  witii  that  subregion  of  the  plate.  By  plotting  loadpath  lines  in  and  around  the 
recirculation  zone,  it  was  possible  to  define  rework  botmdaries  that  essentially 
removed  the  recirculation  zones  at  the  top  and  bottom  of  the  hole. 

After  the  FE  model  of  this  design  was  created  and  analysed,  tiie  loadflow  orientations 
were  then  computed  and  the  results  are  shown  in  Figure  20.  The  amount  of  material 
removed  corresponds  to  a  maximum  rework  depth  of  approximately  2.5  mm  (compare 
the  circular  hole  shape  in  Figme  7  with  the  reworked  hole  shape  in  Figure  20).  The 
removal  of  material  has  resulted  in  a  vertical  shift  in  the  location  of  the  recirculation 
zone,  and  the  loadflow  orientations  are  much  the  same  as  those  obtained  for  the 
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original  circular  geometry.  Also,  because  of  the  localised  nature  of  the  rework  region, 
there  is  no  significant  change  in  the  maximiun  SCF  at  the  edge  of  the  hole  (the  SCF 
increased  by  0.4%).  Hence,  no  benefit  has  been  achieved  at  the  peak  stress  location  by 
the  remov^  of  the  original  recirculation  zone.  The  stresses  in  the  reworked  region  are 
altered  from  their  initial  compressive  stress  state  to  a  state  of  tension,  but  the  stress 
levels  remain  quite  low.  The  loadpath  lines  after  the  removal  of  the  recirculation  zone 
are  shown  in  Figure  21.  Keeping  in  mind  that  an  optimal  rework  shape  has  a  uniform 
state  of  stress  along  most  of  its  bormdary  [1],  the  fact  that  there  is  still  a  considerable 
variation  in  stress  along  the  present  reworked  boimdary  is  indicative  of  the  highly 
non-optimal  nature  of  this  rework  shape. 

4.2  Reworked  circular  hole  in  a  4:1  biaxial  stress  field 

The  loadflow  orientations  obtained  around  a  circular  hole  in  a  large  square  plate 
biaxially  loaded  by  a  4:1  stress  field  have  previously  been  shown  in  Figure  13.  The 
most  significant  recirculation  zone  corresponds  to  the  resolved  ;c-direction  loadflow 
orientations  (Figure  13a)  and  occurs  at  the  top  and  bottom  of  the  hole.  With  the  help  of 
loadpath  lines,  the  rework  region  is  selected  to  coincide  with  that  subregion  of  the 
plate.  A  smaller  recircidation  zone  exists  for  the  resolved  y-direction  loadflow,  and  this 
zone  falls  within  the  a:-direction  recircxilation  zone.  By  plotting  loadpath  lines  in  and 
arotmd  the  recirculation  zone,  it  was  possible  to  define  rework  boxmdaries  that 
essentially  removed  the  zones  of  recircidation  that  occur  at  the  top  and  bottom  of  the 
hole. 

After  the  FE  model  of  this  design  was  created  and  analysed,  the  new  loadflow 
orientations  were  then  computed  and  the  results  are  shown  in  Figure  22.  The  amount 
of  material  removed  corresponds  to  a  maximum  rework  depth  of  approximately  3.0 
mm  (compare  the  circular  hole  shape  in  Figure  13  with  the  reworked  hole  shape  in 
Figure  22).  The  loadpath  lines  are  shown  in  Figure  23.  In  contrast  to  the  xmiaxial 
loading  case,  the  recirculation  region  has  been  cut  out  by  the  rework.  However, 
because  the  amount  of  material  removed  was  fairly  minimal,  as  well  as  being 
restricted  to  the  upper  and  lower  areas  of  the  hole,  there  was  no  significant  change  in 
the  SCF  at  the  edge  of  the  hole  (a  SCF  reduction  of  0.2%).  Although  the  stresses  in  the 
reworked  region  are  altered  from  their  initial  compressive  stress  state  to  a  state  of 
tension,  the  stress  levels  remain  quite  low  compared  to  the  peak  stress  levels  elsewhere 
in  the  structure. 


5.  Discussion  and  Conclusion 

This  report  presents  recent  developments  and  advances  in  the  analysis  and 
visualisation  of  structured  loadflow,  with  the  view  to  establishing  FEA-based 
numerical  procedures  for  loadflow  visualisation  to  gain  insights  into  structural 
performance.  The  particular  eureas  that  have  been  addressed  include:  i)  enhancement  of 
the  theoretical  formulation  and  ntunerical  procedtire,  ii)  studies  of  the  differences 
between  non-optimal  and  optimal  shapes,  and  ili)  the  investigation  of  potential 
applications  to  optimal  shape  reworking. 
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Initially,  the  theory  of  Kelly  and  Elsley  has  been  revisited  and,  in  a  key  development, 
we  have  derived  explicit  equations  for  the  computation  of  loadflow  orientations.  These 
equations  provide  better  accuracy  compared  to  the  previous  approach,  which  was 
based  on  the  iterative  solution  of  non-linear  equations.  Our  derivation  has  also 
dispersed  an  element  of  confusion  arising  from  the  reference  to  the  major  principal 
stress  angle  in  the  presentation  of  the  original  theory.  Then,  by  applying  a  fluidflow 
analogy  proposed  by  Kelly  [16],  we  have  used  the  vector-like  field  of  loadflow 
orientations  to  determine  the  loadpaths  in  a  uniaxiaUy  loaded  square  plate  with  a 
circular  hole.  Here  we  found  that  the  loadflow  orientations  obtained  by  using  the 
explicit  equations  improved  the  fidelity  of  the  calculated  loadpaths.  For  this  case  we 
have  also  studied  the  compaction  of  loadflow  lines  to  determine  the  stress 
concentration  factor,  successfully  verif5dng  the  use  of  the  fluidflow  analogy  in  the 
calculation  of  the  loadpaths.  The  exact  nature  of  this  analogy  is  not  understood  at  this 
point  in  time,  and  it  should  be  studied  in  future  work  to  gain  a  better  understanding. 

The  above  procedures  were  then  applied  to  investigate  the  features  of  non-optimal  and 
optimal  holes  in  plates  where  both  loadflow  orientations  and  loadpaths  were 
determined.  Here  an  optimal  hole  shape  is  taken  to  be  one  that  produces  the  minimum 
peak  SCF  around  the  hole  botmdary.  Through  this  comparison,  it  was  found  that  a 
distinct  recircrdation  zone  is  apparent  for  non-optimal  hole  shapes,  whereas  no 
recirculation  region  is  present  for  optimal  shapes.  In  addition,  the  unportance  of  using 
highly  refined  FE  meshes  to  enable  fine  details  in  recirculation  zones  to  be  resolved 
must  be  emphasised.  Although  very  highly  refined  FE  meshes  were  utilised  in  our 
investigations,  the  implications  are  that  even  more  refined  meshes  are  required  to  fully 
capture  the  complex  behaviour  that  exists  in  recirculation  zones. 

In  the  context  of  shape  optimisation,  the  use  of  loadflow  visualisation  has  been 
investigated  as  a  potential  aid  in  determining  good  initial  estimates  of  rework  shapes. 
Although  the  removal  of  the  recirculation  zone  for  a  non-optimal  hole  shape  leads  to  a 
better  shape,  the  improvement  in  peak  stress  is  insignificant,  even  when  a  distinct 
recirculation  zone  is  present.  This  is  largely  due  to  the  fact  that  the  recirculation  zone  is 
quite  localised  in  nature,  whereas  the  achievement  of  an  optimal  shape  requires  that 
tile  rework  region  be  comprised  of  a  substantial  portion  of  the  original  boundary  in 
order  to  enable  a  significant  reduction  in  SCF  to  be  obtained.  In  situations  where  the 
original  structural  shape  is  similar  to  tiie  optimal  rework  shape,  the  recirculation  zone 
is  less  pronounced  and  it  becomes  more  difficult  to  clearly  identify  a  rework  shape  that 
may  improve  the  design. 

In  general,  the  visualisation  of  loadflow  orientations  and  loadpaths  presents  structural 
designers  and  analysts  with  a  powerful  tool  to  complement  ihose  already  at  their 
disposal.  Loadflow  visualisation  can  improve  their  insight  into  the  transfer  of  load 
through  a  structure,  thus  placing  them  in  a  better  position  to  gauge  how  well  a 
structure  is  performing,  potentially  resulting  in  worthwhile  improvements  in 
structural  efficiency.  The  simplicity  of  our  improved  method  of  calculation  of  loadflow 
orientations  means  that  it  is  now  just  as  easy,  if  not  in  fact  easier,  to  compute  loadflow 
orientations  than  it  is  to  calculate  principal  stresses.  Thus,  there  is  considerable 
potential  to  have  the  computation  of  loadflow  orientations  implemented  as  a  standard 
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feature  in  FEA  codes,  in  much  the  same  way  that  the  calculations  of  principal  stresses 
are  today. 
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Table  1  Stress  concentration  factor  SCF  along  the  x-axis  in  the  vicinity  of  a  circular  hole  of 
radius  a  in  a  uniaxially  loaded  square  plate  computed  from  compaction  of  loadpath 
lines. 


afxc 

Xc  (mm) 

Ac  (mm) 

Xb  (mm) 

Ab  (mm) 

SCF 

0.9524 

10.500 

0.01 

1.3993 

0.0280 

2.8000 

0.9434 

10.600 

0.20 

1.6794 

0.5429 

2.7145 

0.9434 

10.600 

0.01 

1.6819 

0.0269 

2.6900 

0.9346 

10.700 

0.01 

1.9535 

0.0267 

2.6700 

0.9259 

10.800 

0.01 

2.2191 

0.0264 

2.6400 

0.9174 

10.900 

0.01 

2.4779 

0.0258 

2.5800 

0.9091 

11.000 

0.01 

2.7327 

0.0250 

2.5000 

0.8929 

11.200 

0.01 

3.2218 

0.0239 

2.3900 

0.8772 

11.400 

0.01 

3.7065 

0.0238 

2.3800 

0.8624 

11.595 

0.01 

4.1633 

0.0227 

2.2700 

0.8333 

12.000 

0.20 

5.0592 

0.4241 

2.1205 

0.8000 

12.500 

0.20 

6.0912 

0.3998 

1.9990 

0.7692 

13.000 

0.20 

7.0628 

0.3796 

1.8980 

0.7407 

13.500 

0.20 

7.9673 

0.3485 

1.7425 

0.7143 

14.000 

0.20 

8.8200 

0.3283 

1.6415 

0.6667 

15.000 

0.20 

10.4339 

0.3146 

1.5730 

0.6250 

16.000 

0.20 

11.9538 

0.2906 

1.4530 

0.5882 

17.000 

0.20 

13.3776 

0.2838 

1.4190 

0.5556 

18.000 

0.20 

14.7319 

0.2633 

1.3165 

0.5263 

19.000 

0.20 

16.0304 

0.2496 

1.2480 

0.5000 

20.000 

0.20 

17.2718 

0.2496 

1.2480 

0.4545 

22.000 

0.20 

19.7167 

0.2394 

1.1970 

0.4167 

24.000 

0.20 

22.0530 

0.2326 

1.1630 

0.3846 

26.000 

0.20 

24.3433 

0.2257 

1.1285 

0.3571 

28.000 

0.20 

26.5897 

0.2232 

1.1160 

0.3125 

32.000 

0.20 

30.9716 

0.2183 

1.0915 

0.2857 

35.000 

0.20 

34.1956 

0.2155 

1.0775 

0.2778 

36.000 

0.20 

35.2537 

0.2104 

1.0520 

0.2500 

40.000 

0.20 

39.4768 

0.2123 

1.0615 

0.2222 

45.000 

0.20 

44.7117 

0.2102 

1.0510 

0.1961 

51.000 

0.20 

50.9355 

0.2062 

1.0310 

0.1667 

60.000 

0.20 

60.1853 

0.2046 

1.0230 

0.1429 

70.000 

0.20 

70.3166 

0.2002 

1.0010 

0.1250 

80.000 

0.20 

80.3687 

0.2006 

1.0030 

0.1111 

90.000 

0.20 

90.2631 

0.1976 

0.9880 

0.1053 

95.000 

0.20 

95.1443 

0.1976 

0.9880 

0.1005 

99.500 

0.20 

99.6036 

0.1995 

0.9975 
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Table  2  Stress  concentration  factor  SCF  in  the  vicinity  of  a  circular  hole 


of  radius  a  in  a  uniaxially  loaded  large  square  plate  as  computed 
from  finite  element  analysis  nodal  stresses  tahen  along  the  x-axis. 


a/x 

X  (mm) 

SCF 

0.1000 

100.0000 

0.9643 

0.1053 

94.9301 

0.9754 

0.1110 

90.1187 

0.9835 

0.1169 

85.5527 

0.9901 

0.1231 

81.2194 

0.9957 

0.1297 

77.1071 

1.0006 

0.1366 

73.2045 

1.0051 

0.1439 

69.5008 

1.0092 

0.1515 

65.9860 

1.0131 

0.1596 

62.6503 

1.0168 

0.1681 

59.4848 

1.0205 

0.1771 

56.4806 

1.0242 

0.1865 

53.6296 

1.0279 

0.1964 

50.9240 

1.0317 

0.2068 

48.3563 

1.0357 

0.2178 

45.9196 

1.0399 

0.2293 

43.6070 

1.0444 

0.2415 

41.4124 

1.0493 

0.2543 

39.3297 

1.0546 

0.2677 

37.3532 

1.0605 

0.2819 

35.4774 

1.0670 

0.2968 

33.6973 

1.0743 

0.3124 

32.0079 

1.0824 

0.3289 

30.4047 

1.0917 

0.3462 

28.8832 

1.1022 

0.3644 

27.4393 

1.1143 

0.3836 

26.0690 

1.1281 

0.4037 

24.7685 

1.1441 

0.4249 

23.5344 

1.1626 

0.4472 

22.3632 

1.1840 

0.4706 

21.2517 

1.2091 

0.4951 

20.1969 

1.2384 

0.5209 

19.1959 

1.2727 

0.5481 

18.2459 

1.3131 

0.5766 

17.3443 

1.3607 

0.6065 

16.4887 

1.4170 

0.6379 

15.6768 

1.4835 

0.6709 

14.9062 

1.5625 

0.7055 

14.1749 

1.6563 

0.7418 

13.4809 

1.7679 

0.7799 

12.8223 

1.9009 

0.8199 

12.1972 

2.0597 

0.8618 

11.6041 

2.2494 

0.9057 

11.0411 

2.4762 

0.9518 

10.5069 

2.7478 

1.0000 

10.0000 

3.1165 
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Figure  1  Schematic  of  a  force  "stream  tube"  comprised  of  loadpath  walls  proposed  by  Kelly 
and  Elsley  [3]  [4],  shown  for  the  resolved  x-direction  case. 
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(a)  Complete  finite  element  mesh 


(b)  Detail  of  finite  element  mesh  arotmd  circular  hole 


Figure  4  Finite  element  mesh  for  modelling  a  large  square  plate  containing  a  circular  hole. 
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(a)  Orientations  of  principal  stress 


(b)  Orientations  of  principal  stress  a 

Figures  Orientations  of -principal  stresses  a^j  and  around  a  circular  hole  in  a  large 
square  plate  uniaxially  loaded  in  the  y-direction. 
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(b)  Resolved  y-direction  loadflow  orientations 

Figure  7  Resolved  x-direction  and  y-direction  loadflow  orientations  obtained  around  a 
circular  hole  in  a  large  square  plate  uniaxially  loaded  in  the  y-direction. 
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Dimensionless  position  along  x  -axis,  a/x 


Figure  11  Comparison  of  stress  concentration  factors  computed  from  compaction  of  loadpath 
lines  and  finite  element  model  nodal  stresses  for  a  uniaxially  loaded  large  square 
plate  with  a  circular  hole. 
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(b)  Resolved  y-direction  loadflow  orientations 


Figure  13  Resolved  x-direction  and  y-direction  loadflow  orientations  obtained  around  a 
circular  hole  in  a  large  square  plate  under  y:x  =  4:1  biaxial  loading. 
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Figure  14  Resolved  x-direction  and  y-direction  loadpaths  obtained  around  a  circular  hole  in  a 
large  square  plate  under  y:x  =  4:1  biaxial  loading. 
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(a)  Complete  finite  element  mesh 


(b)  Detail  of  finite  element  mesh  aroimd  elliptical  hole 


Figure  17  Finite  element  mesh  for  modelling  a  large  square  plate  containing  a  y:x  =  4:1 
elliptical  hole. 
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(a)  Resolved  a:-direction  loadflow  orientations 


(b)  Resolved  y-direction  loadflow  orientations 


Figure  18  Resolved  x-direction  and  y-direction  loadflow  orientations  obtained  around  an 
optimal  4:1  elliptical  hole  in  a  large  square  plate  under  y:x  =  4:1  biaxial  loading. 
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(a)  Resolved  A:-direction  loadpaths 


Figure  19  Resolved  x-direction  and  y-direction  loadpaths  obtained  around  an  optimal  4:1 
elliptical  hole  in  a  large  square  plate  under  y:x  =  4:1  biaxial  loading. 
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(a)  Resolved  a:-direction  loadflow  orientations 


Figure  20  Resolved  x-direction  and  y-direction  loadflow  orientations  obtained  around  a 
reworked  circular  hole  in  a  large  square  -plate  uniaxially  loaded  in  the  y-direction. 
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(a)  Resolved  x-direction  loadflow  orientations 


Figure  22  Resolved  x-direction  and  y-direction  loadflow  orientations  obtained  around  a 
reworked  circular  hole  in  a  large  square  plate  under  y:x  =  4:1  biaxial  loading. 


41 


DSTO-RR-0166 


Figure  23  Resolved  x-direction  and  y-direction  loadpaths  obtained  around  a  reworked  circular 
hole  in  a  large  square  plate  under  y:x  =  4:1  biaxial  loading. 
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Appendix  A 

Extracting  Stresses  From  NASTRAN  FEA  Results 


A.1  General  comments 

In  order  to  be  able  to  perform  loadflow  calculations  for  a  structure  modelled  in  a  FEA 
package,  it  is  necessary  to  be  able  to  extract  out  the  stresses  that  occur  at  nodes  and 
elements  in  the  FE  model.  These  results  are  saved  in  separate  files  and  used  as  input 
fdes  by  the  program  that  computes  loadflow  orientations.  When  using  the 
MSC.NASTRAN  FEA  package  and  MSC.PATRAN  pre-  and  post-processor,  the 
procedtire  fliat  was  used  to  create  the  nodal  stress  file  is: 

1.  Use  the  Fortran  program  pnode.f  to  pull  node  data  lines  from  the  NASTRAN  data 
deck  and  store  the  results  in  output  file  nds.dat. 

2.  Use  PATRAN  utilities  available  under  the  menu  selection  Utilities/Results/Results 
Utilities  to  create  individual  files  containing  nodal  global  stress  outputs,  one 
component  per  file  (sxx.dat,  syy.dat,  and  sxy.dat). 

3.  Use  the  Fortran  program  snode.f  to  combine  the  above  files  and  calculate  the 
values  of  the  principal  stresses,  and  store  the  results  in  output  file  lfsn.dat. 

A  process  is  used  to  create  the  element  stress  file  containing  the  stresses  at  the 

element  centroids.  The  element  stresses  are  sometimes  used  to  produce  a  less  dense 
display  of  loadflow  orientations  only. 

1.  Use  the  Fortran  program  pelem.f  to  pull  element  lines  from  the  NASTRAN  data 
deck  and  store  the  results  in  output  file  elms.dat. 

2.  Use  PATRAN  utilities  to  create  individual  files  containing  the  element  centroid 
global  stress  outputs,  one  component  per  file  (exx.dat,  eyy.dat,  and  exy.dat). 

3.  Use  flie  Fortran  program  selem.f  to  combine  the  above  files  and  calculate  the 
values  of  the  principal  stresses,  and  store  the  results  in  output  file  to  lfse.dat. 

The  resulting  files  lfsn.dat  and  lfse.dat  are  the  input  files  to  the  program  for  computing 
loadflow  orientations,  which  is  described  in  Appendix  B.  The  different  programs 
referred  to  above  are  listed  in  the  following  four  sections. 

A.2  Fortran  program  pnode.f  for  extracting  nodal  data  lines 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
program  pnode 

C  This  program  writes  all  grid  points  (nodes)  from  the  NASTRAN  data 
C  deck  file  whose  name  is  given  in  the  variable  nameif. 
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c 

C  The  results  are  written  into  the  file  whose  name  is  given  in  the 
C  variable  nameof . 

character  aline*80,bline*80, tokens (5) *30 
character  nameif *80, nameof *80 
integer  index 

nameif  =  'lfOclrwk.dat' 
nameof  =  'lfOclrwknode.dat' 

write(*,*)  'Name  of  input  file  =  nameif (1: index (nameif, '  ')) 
write(*,*)  'Name  of  output  file  =  ', nameof (1: index (nameof,  '  ')) 

open (unit=10 , f ile=nameif ) 
open (unit=20 , f ile=nameof ) 

n  =  0 

do  while  (.true.) 

read (10, ' (a) ' ,end=2000)  aline 
if  (aline (1 :4) .eq. 'GRID')  then 
n  =  n  +  1 

if  (aline (5 :5) .eq. ' * ' )  then 
read (10, ' (a) ' )  bline 
aline=aline (1 : 71) //bline (9 : ) 
call  gettokens (aline, '  ' ,ntokens, tokens) 
call  chkf loat ( tokens ( 3 ] ) 
call  chkfloat (tokens (4) ) 
call  chkfloat (tokens (5) ) 

aline  =  tokens (2) (l:index( tokens (2) , '  '))// 

Sc  tokens(3)  (l:index(tokens(3) ,  '  '))// 

&  tokensU)  (l:index(tokens(4)  ,  '  '))// 

Sc  tokens  (5)  (1:  index  (tokens  (5) ,  '  ')) 

read(aline, *)  nid,x,y, z 
else 

read (aline, ' (8x, i8, 8x,3f8 .4) ' )  nid,x,y, z 
endif 

write(20, ' (a, iS, a, 3 (a, fl2 . 4) ) ') 

Sc  'GRID,  ',nid,  ,x,  ','  ,y,  ’,  ',z 

end  if 
end  do 

2000  continue 

write (*,*)  'Number  of  GRID  points  processed  =  ',n 

close (10) 
close (20) 

stop 

end 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCGGCCCCCC 

subroutine  chkfloat (s) 

character  s*(*) 

integer  index 

n  =  index (s, '  ' ) -1 
i  =  n-2 

if  (i.gt.l.and. (s(i:i) .eq. '-'  .or.  s (i: i) .eq. ' +' ) )  then 
s  =  s(l:i-l)//'E'//s(i:) 
endif 

return 

end 

GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCGGGGGGGGG 
subroutine  gettokens (string, delimiters, ntokens, tokens) 
implicit  none 
integer  ntokens 

character  string* (*), delimiters* (*), tokens (*)* (*) 
character  token* 80 
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character  strtok*80,char*l 
ntokens  =  0 

tokens (1)  =  strtok (string,  delimiters) 

do  while  (token  .ne.  char(O)) 
ntokens  =  ntokens  +  1 
token  =  strtok(char (0) ,  delimiters) 
if  (token.ne.char (0) )  then 
tokens (ntokens+l)  =  token 
endif 
enddo 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

character*80  function  strtok (srce,  delim) 

Tokenize  a  string  in  a  similar  manner  to  C.  The  usage  is 
very  similar  to  the  C  function. 

Input:  srce  =  source  string  to  tokenize  (see  usage  note) . 
delim  =  delimiter  string,  used  to  determine  the 
beginning/end  of  each  token  in  a  string. 

Output:  strtok  =  string  token. 

Usage:  First  call  strtok  with  the  string  to  tokenize 

as  srce,  and  the  delimiter  string  used  to  tokenize. 

Then  call  strtok  with  srce  set  to  ohar(O) . 

If  strtok  returns  char(O) ,  no  more  tokens  are 
available . 

implicit  none 

character  srce*(*),  delira* (*) 

character  savestr*1024 

integer  start.  In,  begpos,  endpos 

save  start,  savestr 

integer  index 
character  char*l 

if  (srce(l:l)  .ne.  char(O))  then 
start  =  1 

if  (len(srce) .gt.len(savestr) )  then 

write (*,*)  '***  warning:  source  string  too  long  for  buffer.' 
endif 

savestr  =  srce 
endif 

begpos  =  start 
In  =  len (savestr) 

5  continue 

if  ((begpos  .le.  In)  .and. 

&  (index (delim, savestr (begpos :begpos) )  .ne.  0)  )  then 
begpos  =  begpos  +  1 
goto  5 
endif 


if  (begpos  .gt.  In)  then 
strtok  =  char(O) 
return 
endif 

endpos  =  begpos 

10  continue 

if  ((endpos  .le.  In)  .and. 

Si  (index  (delim,  savestr  (endpos  :endpos) )  .eq.  0)  )  then 
endpos  =  endpos  +  1 
goto  10 
endif 

strtok  =  savestr(begpos:endpos-l) 
start  =  endpos 
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return 

end 


A.3  Fortran  program  pelem.f  for  extracting  element  data  lines 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

program  pelem 

C  This  program  writes  all  elements  from  the  NASTRAN  data  deck  file 
C  whose  name  is  given  in  the  variable  nameif. 

C 

C  The  results  are  written  into  the  file  whose  name  is  given  in  the 
C  variable  nameof. 

character  al ine*  8  0 , nameif *  8  0 , nameof  *  8  0 
integer  index 

nameif  =  'lfOclrwk.dat' 
nameof  =  'lfOclrwkelem.dat' 

write(*,*)  'Name  of  input  file  =  nameif (1 : index(nameif,  '  ')) 
write (*,*)  'Name  of  output  file  =  ' .nameof (1: index (nameof, '  ')) 

open (unit=10, f ile=nameif ) 
open (unit=20 , f ile=nameof ) 

n  =  0 

do  while  ( . true . ) 

readdO,  '  (a)  '  ,end=2000)  aline 
if  (aline (1:6) .eg. 'CQUAD4 ' )  then 
n  =  n  +  1 
i  =  len (aline) 

do  while  (aline(i:i) .eg. '  ’) 
i  =  i  -  1 
end  do 

write  (20,  '  (a)  ' )  aline  (l:i) 
end  if 
end  do 

2000  continue 

write(*,*)  'Number  of  elements  processed  =  ',n 

close (10) 
close (20) 

stop 

end 


A.4  Fortran  program  snode.f  for  combining  nodal  stresses 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

program  snode 
character  aline*80 

open (unit=10 , f ile= ' If Oclrwknode . dat ' ) 
open(unit=20,file='lf0olrwknxx.dat' ) 
open (unit=30 , f ile= ' lf0clrwknyy.dat ' ) 
open(unit=40, file=' lf0clrwknxy.dat ' ) 
open{unit=50, file= ' lfOclrwksn.dat' ) 

do  i=l,18 
read(20,*) 
read(30,*) 
read  (40,  *) 
end  do 

write (50, ' (a8, 2x, 5 (al0,2x) ) ' ) 

&  ' id' , 'x' , 'y' , ' sxx’ , ' syy' , ' sxy' 
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n  =  0 

do  while  ( . true . ) 

readdo,  '  (a)  '  ,eiid=10)  aline 
aline  =  aline (6: 13)//'  ■ //aline (16 : ) 
read(aline, *)  id,x,y,z 
read(20,*>  iduntmy,  j dummy,  sxx 
read(30,*)  idummy, j dummy, syy 
read{40,*)  idummy, j dummy, sxy 


n  =  n+1 

write(S0, ' (i8,5(2x,fl0.4)) ')  id, x,y, sxx, syy, sxy 
end  do 

10  write(*,*)  'Number  of  nodes  processed  =  ',n 

close (10) 
close (20) 
close (30) 
close  (40) 
close (50) 

stop 

end 


A.5  Fortran  program  selem.£  for  combining  element  centroid  stresses 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
program  selem 
character  slab* 8 

open (unit=l0 , f ile= '  lf0clrwkelem.dat ' ) 
open  (unit=20 ,  f  ile= '  If  0clrw)texx .  dat ' ) 
open  (unit=30 ,  f  ile= '  If  Oclrwjceyy .  dat  ’ ) 
open (unit=40, file=' If 0clrwkexy.dat ' ) 
open {unit=50 , f ile= ' If Oclrwkse . dat ' ) 

do  i=l,18 
read(20,*) 
read(30,  *) 
read(40,  *) 
end  do 

write (50, ' (5a8,3 (2x,al0) ) ' )  'id' , 'nl' , 'n2' , 'n3' , 'n4' , 
t  'sxx' , 'syy' , 'sxy' 

n  =  0 

do  while  (.true.) 

read ( 10 , * , end=l 0 )  slab , id , n j unk , nl , n2 , n3 , n4 
read(20,*)  idummy, j dummy, sxx 
read(30,*)  idummy, j dummy, syy 
read (40 , * )  idummy, jdummy, sxy 

n  =  n+1 

write(50, ' (5i8 , 3 (2x, flO .4) ) ')  id,nl,n2 ,n3 ,n4 , sxx, syy, sxy 
end  do 

10  write(*,*)  'Number  of  elements  processed  =  ',n 

close (10) 
close (20) 
close (30) 
close (40) 
close (50) 

stop 

end 
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Appendix  B 

Determining  Loadflow  Orientations 


B.l  General  comments 

A  Fortran  program  lfvec±or.f  was  developed  for  determining  the  loadflow  orientations 
from  a  set  of  datafiles  containing  node  and  element  stresses.  In  our  case,  these 
intermediate  datafiles  are  obtained  by  using  the  MSC.PATRAN  FEA  post-processing 
software  to  extract  the  desired  stresses  from  die  FE  analysis  performed  by 
MSC.NASTRAN.  The  results  produced  by  lfvector.f  are  written  to  an  output  file  in  a 
Tecplot-compatible  format  that  enables  Tecplot  to  be  used  for  visualisation  of  the 
loadflow  orientations  for  the  structure  that  has  been  analysed.  Both  the  x-direction  and 
y-direction  loadflow  orientations  are  computed,  and  the  results  are  stored  in  two 
formats:  as  a  field  of  unit  vectors  with  x-  and  y-components  set  to  produce  the  required 
angular  orientation  (for  use  within  Tecplot),  and  as  angles  in  degrees.  The  program 
also  computes  the  principal  stress  directions  and  stores  them  in  a  similar  manner, 
enabling  comparisons  widi  loadflow  orientations  to  be  easily  obtained.  A  hsttng  of  the 
lfvector.f  program  is  provided  in  the  next  section. 


B.2  Fortran  program  Ifvectori  for  calculating  loadflow  orientations 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
program  Ifvector 
c  10t.06-1998  Witold  Waldman 

c 

C  Th.is  program  computes  the  orientations  of  loadpath.  vectors  from  a 
C  set  of  global  stress  data  (sjdc,  syy,  and  sxy)  at  a  series  of  points. 

Q 

C  It  then  outputs  the  results  of  the  calculations  in  Tecplot  format. 

C 

C  17-08-1998  Witold  Waldman 

C  ,  ,  , 

C  Modified  code  to  accept  file  names  from  coiranand  line.  This  is 
C  accomplished  by  making  use  of  HP-UX  operating  system  intrinsic 
C  functions  and  subroutines  that  are  available  in  the  Fortran  language. 

C  To  link  with  the  required  support  libraries,  it  is  necessary  to 
C  supply  the  +U77  command  line  switch  to  the  Hewlett-Packard  f77 
C  compiler. 

C 

C  21-05-1999  Witold  Waldman 

Q 

C  Modified  code  by  updating  the  equations  used  to  compute  loadflow 
C  orientations  to  the  latest  theoretical  formulation. 

implicit  none 

integer  lui,  luo,maxnode,maxelem,maxnptr 

parameter  (lui=l) 
parameter  (luo=2) 
parameter  (maxnode=100000) 
parameter  (maxelem=100000) 
parameter  (maxnptr=100000) 

character  nodefile* 80 
character  elemfile*80 
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character  vectorfile*80 
integer  i,j 

integer  ios 

integer  nelem 

integer  nnode 

integer  nlist 

integer  nclarg 

integer  lenn 

integer  lene 

integer  lenv 

integer  nodeptr (maxnptr) 

integer  nodenum(maxnode) 

integer  nodelist (maxnode) 

real*8  nodex (maxnode) 

real*8  nodey (maxnode) 

real*8  nodesxx (maxnode) 

real*8  node syy (maxnode) 

real*8  nodesxy (maxnode) 

real*8  nodesll (maxnode) 

real*8  nodes22 (maxnode) 

real*8  nodetll (maxnode) 

real *8  nodet22 (maxnode) 

real *  8  nodexsn (maxnode ) 

realms  nodextau (maxnode) 

real*8  nodexalpha (maxnode) 

real *  8  nodeysn (maxnode ) 

real*8  nodeytan (maxnode) 

real *  8  nodeyalpha (maxnode ) 

real*8  listsll (maxnode) 

real*8  lists22 (maxnode) 

real*8  listtll (maxnode) 

real*8  listt22 (maxnode) 

real*  8  listxsn(mcoaiode) 

real*8  listxtau (maxnode) 

real*8  listxalpha (maxnode) 

real *8  listysn (maxnode) 

real*8  listytau (maxnode) 

real*8  listyalpha (maxnode) 

integer  elemnum (maxelem) 

integer  elemnl (maxelem) 

integer  elemn2  (mcixelem) 

integer  elemnS (maxelem) 

integer  elenin4  (maxelem) 

real*8  elemxc (maxelem) 

real*8  elemyc (maxelem) 

real*8  elemsxx (maxelem) 

real*8  elerasyy (maxelem) 

real*8  elemsxy (maxelem) 

real*8  elemsll (maxelem) 

real*8  elems22 (maxelem) 

real*8  elemtll (maxelem) 

real*8  elemt22 (maxelem) 

real*8  elemxsn  (mcocelem) 

real*8  elemxtau  (meutelem) 

real*8  elemxalpha (maxelem) 

real*8  elemysn (maxelem) 

real*8  elemytau (maxelem) 

real*8  elemyalpha (maxelem) 

real*8  xl,x2,x3,x4 

real*8  yl,y2,y3,y4 

logical  doselectednodes 

logical  doelements 

logical  donodes 

real *8  cosd 

real* 8  sind 

integer  jmod 

integer  iargc 

!  Specify  the  names  of  the  data  files  that  contain  stress  results 
!  for  individual  nodes  and  elements,  as  well  as  the  name  of  the 
!  output  file. 

nolarg= iargc ( ) 

if  (nclarg. eg. 3)  then 
call  getarg (1, nodefile) 
call  getarg (2, elemfile) 
call  getarg(3,vectorfile) 
else 
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write(*, '  (/,lx,a,/) ')  _ 

&  'Usage:  Ifvector  nodalstressfile  elementstressfile  outputfile' 
stop 
endif 

call  lenstr(iiodefile,lemi) 
call  letLStr(elemfile,lene) 
call  leiistr(vectorfile,lenv) 


v7rite(*,*) 

writet*,*)  'Pile  of  node  stresses  (input)  :  ' ,nodefile(l:lenn) 

writei*,*)  'File  of  elem  stresses  (input)  :  ' ,elemfile(l:lene) 

write(*,*)  'File  of  load  vectors  (output):  ' ,vectorfile(l:lenv) 

write  {*,’‘) 

doselectednodes= . false . 
doelement  s= . true . 
donodes=.true. 

!  Read  in  the  stress  results  for  the  centre  of  each  element. 

open (unit=lui, f ile=elemfile, status= ' old' ) 

i=0 

ios=0 

readdui,  *) 
do  while  (ios.eq.O) 
i=i+l 

read(lui,*,iostat=ios)  elenmum(i) , 

elemnl  (i)  ,eleTnn2  (i)  ,elemn3  (i)  ,elemn4  (i) , 
elemsxx(i)  ,elemsyy(i)  ,elemsxy(i) 

if  (ios.eq.O)  then 

call  sls2tlt2 (elemsxx(i) ,elemsyy(i) ,elemsxy (i) , 

elemsll(i) ,elems22 (i) ,elemtll(i) ,elemt22 (i) ) 

endif 
enddo 
close (lui) 
nelem=i-l 

write(*,’^)  'Number  of  elements  =  '.nelem 
!  Read  in  the  stress  results  at  the  nodes. 


open (unit=lui, file=nodef ile, status= ' old' ) 
i=0 


ios=0 

readdui,*) 
do  while  (ios.eq.O) 
i=i+l 

read (lui, *, iostat=ios) 


nodenum ( i ) , nodex ( i ) , nodey ( i ) , 
nodesxx(i)  ,nodesyy(i)  ,nodesxy(i) 


if  (ios.eq.O)  then 

nodeptr (nodenum (i) ) =i 

call  sls2tlt2 (nodesxx(i) ,nodesyy (i) ,nodesxy (i) , 

nodesll (i) ,nodes22 (i) ,nodetll (i) ,nodet22 (i) ) 


endif 


enddo 
close (lui) 
nnode=i-l 

write (*,*)  'Number  of  nodes  =  ',nnode 


!  Compute  and  store  the  coordinates  of  the  centre  of  each  element. 


do  i=l, nelem 

xl=nodex (nodeptr (elemnl ( i) ) ) 
x2=nodex (nodeptr (elemn2 (i) ) ) 
x3=nodex (nodeptr (elemnS (i) ) ) 
x4=nodex (nodeptr (elemn4 (i) ) ) 
yl=nodey (nodeptr (elemnl (i) ) ) 
y2=nodey (nodeptr (elemn2 (i) ) ) 
y3=nodey (nodeptr (elemnl ( i) ) ) 
y4 =nodey (nodeptr (elemn4 (i) ) ) 
elemxc (i) = (xl+x2+x3+x4) /4 . OdO 
elemyc (i) = (yl+y2+y3+y4) /4 . OdO 
enddo 


!  Compute  angle  of  the  loadpath  vector  for  selected  nodes . 
!  This  feature  may  be  of  assistance  when  debugging. 


if  (doseleotednodes)  then 
nlist  =  2 
nodelist (1)  =  7 
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nodelist (2)  =  1 

call  wrtaxismsgCO,  '  for  selected  nodes') 
do  i=l,nlist 

j=nodeptr (nodelist (i) ) 

write ( * , * )  ' Processing  node ' , nodenum ( j ) ,  ' . ' 
call  sls2tlt2  {nodesxx(j)  ,nodesyy(j)  ,nodesxy(j)  , 

&  listsll (i) , lists22 (i) , listtll (i) ,  listt22 (i) ) 

call  sntau(l,nodesxx(j) .nodesyyCj) ,nodesxy(j), 

&  listxsn(i) ,listxtau(i) ,iistxalpha(i) ) 

call  sntau(2,nodesxx(j)  ,nodesyy(j)  ,nodesxy(j) , 

&  listysn(i)  ,listytau(i)  ,listyalpha(i) ) 

enddo 
endif 

!  Compute  angle  of  the  loadpath  vector  at  each  of  the  nodes, 
if  (donodes)  then 

call  wrtaxismsg(0, '  for  all  nodes') 
do  i=l,nnode 

call  sntau(l,nodesxx(i) ,nodesyy(i) ,nodesxy(i) , 

Sc  nodexsn  ( i ) ,  nodext  au  ( i ) ,  nodexalpha  ( i ) ) 

call  sntau(2,nodesxx(i) ,nodesyy(i) ,nodesxy(i) , 

Sc  nodey sn  ( i ) ,  nodeyt  au  ( i ) ,  nodeyalpha  ( i ) ) 

if  (jmod(i, 10000) .eg. 0)  write(*,*)  i, '  nodes  processed. ' 
enddo 

if  ( jmod(nnode, 10000) .ne.O)  write(*,*)  nnode, '  nodes  processed. ' 
endif 

!  Compute  angle  of  the  loadpath  vector  at  the  centre  of 
!  each  of  the  elements . 

if  (doelements)  then 

call  wrtaxismsg(0, '  for  all  elements') 
do  i=l,nelem 

call  sntau(l,elemsxx(i)  ,elemsyy(i)  ,elemsxy(i) , 

Sc  eiemxsn(i)  ,eiemxtau(i)  , elemxalpha (i)  ) 

call  sntau(2,elemsxx(i) ,elemsyy(i) ,elemsxy(i) , 

Sc  elemysn(i)  ,elemytau(i)  , elemyalpha (i)  ) 

if  { jmod (i, 10000) . eq. 0)  write(*,*)  i,  '  elements  processed.' 
enddo 

if  ( jmod (nelem, 10000) .ne.O)  write(*,*)  nelem. 

Sc  '  elements  processed. ' 

endif 

!  Write  the  results  of  the  stress  calculations  to  an  output 
!  file  in  Tecplot  format. 

if  (donodes. or. doelements. or. doselectednodes)  then 
write (*,*) 

Sc  'Writing  loadflow  orientations  to  Tecplot  output  file.' 
open  (\init=luo,  f  ile:=vectorfile) 
writeduo,  '  (a) ') 

Sc  'titlec="Loadflow  orientations  and  stress  results'" 
write (luo, ' (a) ' )  'variables=' 
write  duo, ' {a6,l9al3) ') 

' "id" ' , ' "X" ' , ' "y" ' , ' "sxx" ' , ' "syy" ' , ' "sxy"  ' , 

' "sll" ' , ' "thetall" ' , ' "sllx" ' , ' "slly" ' , 

' "s22" ' , ' "theta22" ' , ' "s22x" ' , ' "s22y" ' , 

' "xlfth" ' , ' "xlfthx" ' , ' "xlfthy" ' , 

' "ylfth" ' , ' "ylfthx" ' , ' "ylfthy" ' 
if  (doselectednodes)  then 

write ( * , * )  ' Zone :  selected  nodes . ' 
write  duo, ' (a) ' )  'zone  t="Seleoted  nodes"' 
do  i=l,nlist 

j=nodeptr (nodelist (i) ) 
writeduo,  '  (i6,19fl3.4)  ') 

nodenum ( j ) , nodex ( j ) , nodey ( j ) , 
nodesxx ( j ) , nodesyy ( j ) , nodesxy ( j ) , 

listsll (i) , listtll (i) , cosd (listtll (i) ) , sind (listtll (i) ) , 
lists22  (i) ,  listt22  (i) ,  cosd(listt22  (i) ) ,  sinddistt22  (i) )  , 
listxalpha(i) ,cosd(listxalpha(i) ) , sind(listxalpha (i) ) , 
listyalpha(i) ,oosd(listyalpha(i) ) , sind (listyalpha (i) ) 
enddo 
endif 

if  (donodes)  then 

write (*,*)  'Zone:  all  nodes.' 
write  duo, ' (a, i6,a, ie,a) ' ) 

'zone  t="Nodes"  N=', nnode, '  E=', nelem, 

'  FcccFEPOINT  ET=QUADRIIiATERAL' 
do  i=l, nnode 
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& 

& 

& 

& 

£c 

& 


& 


& 

& 

Sc 

& 

& 

& 


writeCluo, ' (i6,19fl3.4)  ' ) 

iiodenum(i)  ,iiodex(i)  ,iiodey{i) , 
nodesxx(i) ,nodesyy(i) ,nodesxy(i) , 

nodeslKi)  .nodetll  (i) ,  cosd(nodetll  (i) )  ,sind(iiodetll{i) ) , 
iiodes22  (i)  ,nodet22  (i) ,  cosd(nodet22  (i) )  ,sind(iiodet22  (i) ) , 
nodexalpha  ( i ) ,  cosd  (nodexalpha  ( i ) ) ,  sind  (nodexalpha  ( i ) ) , 
nodeyalphaii)  ,cosd(nodeyalpha(i) )  ,sind(nodeyalpha{i) ) 
enddo 


do  i=l,nelem 

write (luo,  '  (4i8)  ') 


nodeptr  ( eleitinl  ( i ) )  ,  nodeptr  ( eleran2  ( i ) )  , 
nodeptr  (eletnnS  (i) ) , nodeptr  (elemn4  (i) ) 


enddo 


endif 

if  (doelements)  then 

write ( * , * )  ' Zone :  all  elements . ' 

write (luo, ' (a) ' )  ' zone  t="Elements"  ' 
do  i=l,nelem 

writeduo,  '  (i6,19fl3.4)  ') 
i,elemxc(i) ,elemyc(i) , 
elemsxx(i)  ,elemsyy(i)  ,elemsxy(i) , 

elemsll(i) ,elemtll(i) , oosd(elemtll (i) )  ,  sind(elemtll (i) ) , 
elems22{i) ,elemt22(i) , cosd(elemt22 (i) )  ,  sind(elemt22 (i) ) , 
elemxalpha(i)  ,cosd(elemxalpha(i) ) ,  sind(elemxalpha(i) )  , 
elemyalpha  (i) ,  cosd  (elemyalpha  (i) ) ,  sind  (elemyalpha  (i) ) 
enddo 


endif 
close  (luo) 
endif 


stop 

.end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

subroutine  wrtaxisrasg(iaxis,msg) 


implicit  none 

integer  iaxis 
character  msg*(*) 

if  (iaxis. eg. 1)  then  , 

write ( * , * )  ' Computing  x-direction  loadpath  orientation  / / 

&  msg// '  . ' 

else  if  (iaxis. eq. 2)  then  . 

write (*,*)  'Computing  y-direction  loadpath  orientation'// 
&  msg// ' . ' 

6XS€ 

write (*,*)  'Computing  x-  Sc  y-direction  loadpath  '// 

&  ' orientation' //msg/ / ' . ' 

endif 


return 

end 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

subroutine  sls2tlt2 (sxx, syy, sxy,  si , s2 , thetal , theta2 ) 

C  This  subprogram  computes  the  principal  stresses  si  and  s2  oriented 
C  at  angles  thetal  and  theta2  using  the  stresses  sxx,  syy,  and  sxy. 

c 

C  si  is  the  numerically  largest  principal  stress  and  thetal  is  its 
C  associated  cuigle. 

C  •  .  ,  .  . 

C  s2  is  the  numerically  smallest  principal  stress  and  theta2  is  its 

C  associated  2uigle. 

C 

C  Note  that,  by  the  definition  of  the  principal  angles  obtained  from 
C  "Mechanics  of  Materials"  by  Higdon  et  al.,  one  of  the  principal 
C  stress  angles  lies  between  +/-45  degrees  (inclusive) ,  and  the  other 
C  is  90  degrees  greater. 

implicit  none 

real *8  sxx, syy, sxy, si, s2, thetal, theta2 
real’‘8  psl,ps2,pthl,pth2 
real*8  datand,dcosd,dsind 
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if  (sxx-syy.ne.O.OdO)  then 
pthl  =  datan<i(2  .  OdO*sxy/ (sxx-syy) ) /2  .  OdO 
else 

if  (sxy.gt.O.OdO)  then 
pthl  =  +4 5. OdO 
else  if  (sxy.lt .0 . OdO)  then 
pthl  =  -45. OdO 
else 

pthl  =  O.OdO 
endif 
endif 

pth2  =  pthl  +  9 O.OdO 

if  (pth2.gt.90.0d0)  pth2  =  pth2  -  180. OdO 

psl  =  (sxx+syy) /2 . OdO  +  (sxx-syy) /2 .OdO*dcosd (2 . OdO*pthl)  + 

&  sxy*dsind(2 . OdO*pthl) 

ps2  =  (sxx+syy) /2 . OdO  +  (sxx-syy) /2 . OdO *dcosd (2 . OdO *pth2)  + 

&  sxy*dsind (2 . 0d0*pth2) 

if  (psl.ge.ps2)  then 
si  =  psl 
s2  =  ps2 
thetal  =  pthl 
theta2  =  pth2 
else 

si  =  ps2 
s2  =  psl 
thetal  =  pth2 
theta2  =  pthl 
endif 

return 

end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

subrout ine  snnsnt(sxx,syy,sxy, theta, snn,snt) 

C  This  subprogram  computes  the  normal  eind  shearing  stresses  snn  Md  snt 
C  on  an  arbitrary  plane  through  a  point  oriented  at  an  angle  theta 
C  with  respect  to  a  reference  x  axis  and  the  Jcnown  stresses  sxx,  syy 
C  and  sxy. 

implicit  none 

real *  8  sxx , syy , sxy , theta , snn , snt 

real*8  sintheta, costheta, sin2theta, cos2theta 

real *  8  ds ind , dcosd 

sintheta  =  dsind (theta) 
costheta  =  dcosd (theta) 
sin2theta  =  sintheta**2 
cos2theta  =  costheta**2 

snn  =  sxx*cos2theta  +  syy*sin2theta  +  2 . OdO*sxy*sintheta*oostheta 

snt  =  - (sxx-syy) *sintheta*oostheta  +  sxy* (oos2theta-sin2theta) 

return 

end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

subroutine  sntau ( iaxis , sxx, syy, sxy, sn, tau, alpha) 

C  This  subprogram  determines  the  angle  alpha  which  causes  the  x-  or 
C  y-components  of  the  loadflow  to  be  zero.  The  associated  values  of 
C  normal  stress  sn  cuid  shear  stress  tau  are  also  computed  for  the 
C  estimated  value  of  alpha. 

implicit  none 

integer  iaxis 

real *  8  sxx , syy , sxy , sn , tau , alpha 
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real*8  theta 
real*8  dat2ui2d 

if  (iaxis.ne.l  .and.  iaxis.ne.2)  then  , 

write(*,*)  '***  Invalid  value  for  iaxis  was  supplied.' 
stop 
endif 

if  ( iaxis. eq.l)  then 

!  Consider  force  equilibrium  in  x-direction. 
alpha  =  datan2d(sxy, sxx) 
else  if  (iaxis. eq. 2)  then 

!  Consider  force  equilibrium  in  y-direction. 
alpha  =  datan2d(syy, sxy) 
endif 

theta  =  alpha  -  90.0d0 

call  snnsnt (sxx, syy, sxy, theta, sn, tau) 

return 

end 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

subroutine  lenstr(s,ls) 

C  Determine  the  length  Is  of  a  string  s. 

implicit  none 

character  s*(*) 
integer  Is 

integer  i 

!  Determine  the  location  of  the  rightmost  non-blank  character 
!  and  use  this  to  define  the  length  of  the  string. 

Is  =  len(s) 

10  if  (s(ls:ls) .eq. '  ')  then 
Is  =  Is  -  1 
if  (Is.gt.O)  goto  10 
endif 

!  Check  to  see  if  a  null  character  is  present,  and  adjust  the 
!  length  of  the  string  to  exclude  this. 

do  i=l,ls 

if  (ichar(s(i:i) ) .eq.O)  then 
Is  =  i  -  1 
goto  20 
endif 
enddo 

20  continue 


!  String. 

!  Number  of  characters  in  string. 


return 


Appendix  C 

Loadpath  Visualisation  Using  Tecplot 


The  following  procedure  has  been  successfully  used  for  loadpath  visualisation  using  the 

Tecplot  data  visualisation  program  [18] : 

1.  Conduct  a  FEA  of  the  structure  of  interest  and  compute  the  stresses  a^,  cr^  and  at 
each  of  the  nodes  in  the  FE  model.  Store  the  nodal  coordinates  and  corresponding  nodal 
stresses  in  a  file  for  subsequent  post-processing. 

2.  Use  the  previously  computed  stresses  to  resolve  the  loadflow  orientations  in  the  x-  and 
y-directions  by  using  Equation  (11)  and  Equation  (12).  The  loadflow  orientation  can  be 
represented  by  unit  vectors  with  appropriate  x-  and  y-components  (direction  cosines). 

3.  Produce  a  data  fUe  in  Tecplot  FEPOINT  format  (see  [18],  Chapter  16,  Working  with 
Finite-Element  Data)  where  the  {x,  y)  coordinates  of  the  nodal  points  are  stored,  together 
with  the  direction  cosines  associated  widi  the  computed  loadflow  orientation  at  each  of 
the  points.  This  defines  a  vector-like  field  of  tmit  vectors  with  their  associated 
orientations. 

4.  Open  Tecplot  and  read  in  direction  cosines  of  the  loadflow  orientation,  and  then  activate 
the  Vector  layer  in  the  Tecplot  sidebar. 

5.  Set  the  a:-  and  y-components  of  the  vector  field.  This  is  accomplished  by  choosing  Vector 
Variables  from  the  Field  menu. 

6.  Once  the  Vector  checkbox  has  been  selected  and  the  vector  components  have  been 
chosen,  you  can  redraw  the  Tecplot  screen  and  see  the  vector  plot. 

7.  If  desired,  it  is  possible  to  alter  vector  plot  attributes  using  the  Vector  Attributes  dialog. 
Tecplot  allows  four  different  types  of  vectors:  tail  at  point,  anchored  at  midpoint,  head 
at  point,  and  head  only.  Anchored  at  midpoint  is  a  good  choice.  Now,  because  loadflow 
orientations  are  not  a  true  vector  quantity,  it  is  useful  to  remove  the  arrowhead  from  the 
vector  display.  This  can  be  accomplished  by  setting  the  size  of  the  vector  arrowheads  to 
zero  by  choosing  Vector  Arrowheads  from  the  Field  menu. 

8.  In  Tecplot,  it  is  also  possible  to  modify  the  characteristics  of  many  other  aspects  of  vector 
plots  (e.g.  vector  length).  For  further  details  see  [18],  Chapter  11,  Creating  Vector  Plots. 

9.  Having  defined  the  vector  components  to  be  associated  with  the  field  of  loadflow 
orientations,  you  can  place  streamtraces  one  at  a  time  or  in  groups  called  rakes  (see  [18], 
Chapter  12,  Streamtraces).  A  streamtrace  rake  is  a  set  of  streamtraces  with  starting 
points  along  a  given  line.  You  can  place  streamtraces  or  streamtrace  rakes  by  choosing 
Streamtrace  Placement  from  the  Field  menu.  This  brings  up  the  Streamtrace  Placement 
dialog.  To  place  a  rake  of  streamlines.  Select  the  Enter  Rake  Positions  checkbox  and 
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enter  the  end  positions  for  the  rake  as  {x,  y)  coordinates.  After  defining  where 
streamtraces  are  to  start,  press  the  Place  Streamtrace(s)  button  to  have  them  drawn  on 
the  plot.  If  desired,  you  toggle  the  display  of  the  loadflow  (vector)  orientations  by 
clicking  the  Vector  button  in  the  Tecplot  sidebar. 

10.  In  some  cases,  streamtraces  will  spiral  rather  than  producing  the  desired  closed  contoru*. 
It  is  possible  to  use  a  streamtrace  termination  line  that  terminates  any  streamtraces  that 
cross  it.  To  create  a  streamtrace  termination  line,  from  the  Term  Line  page  of  the 
Streamtrace  Details  dialog,  click  the  Draw  Stream  Term  Line  button.  You  can  then  use 
the  mouse  to  click  the  desired  starting  point  for  the  termination  hne.  Then  click  at 
addition  points  to  define  the  desired  polyline.  To  end  the  termination  Une,  press  the  Esc 
key.  Note  that  only  one  termination  line  can  exist  at  any  one  time  in  a  given  frame.  You 
can  control  the  streamtrace  termination  line  from  die  Term  Line  page  of  the  Streamtrace 
Details  dialog  (e.g.  activating  and  deactivating). 

Tecplot,  uses  a  predictor-corrector  integration  algorithm  to  calculate  streamtraces  [18].  The 
basic  idea  is  to  create  the  streamtrace  by  moving  in  a  series  of  small  steps  from  the  starting 
point  in  the  direction  of,  or  in  opposition  to,  the  local  vector  field.  Each  step  is  only  a 
fraction  of  a  cell  or  element.  Tecplot  automatically  adjusts  the  step  size  based  on  local  cell 
shape  and  vector  field  variation.  The  streamtrace  integration  procedxure  can  be  controlled  by 
modifying  parameters  in  the  Integration  page  of  the  Streamtrace  Details  dialog. 
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